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ABSTRACT 


Bevelopment of quantitative analytical procedures for relating 
the water quality parameter to the characteristics of the back- 
scattered signals, measured by a remote sensor, necessitates further 
rphysical insight in the area of radiative transfer processes in 
turbid media. The present report discusses the applications of a 
Honte 'Carlo simulation .model for radiative transfer in turbid water. 

The model is designed to calculate the characteristics of the back- 
scattered signal from an illuminated body of water as a function of 
the turbidity level, and the spectral properties of the suspended 
particulates. The «optical properties of the environmental waters, 
necessary for model applications., have been derived from available 
experimental data and/or calculated from Mie formalism. Results of 
applications of the model, which have been implemented in support 
of a laboratory program at MSA/'Langley Research Center, are presented. 



TABLE OF CONTENTS 


Page 

LIST OF ILLUSTRATIONS vii 

LIST OF TABLES ix 

1.0 INTRODUCTION 1 

1.1 Optical Parameters of Turbid Water 5 

1.2 Radiative Transfer Model 5 

1.2.1 Monte Carlo Simulation for Narrow Beam Transmission 6 

1.3 Conclusions and Organization of This Report 7 

2.0 MEASURED SCATTERING FUNCTIONS 11 

2.1 Variation of Scattering Function with Wavelength \ 14 

2.2 Scattering Probability Function F(0) 16 

3.0 CALCULATED SCATTERING FUNCTIONS 25 

3.1 Mie Theory for Single Particle Scattering 26 

3.2 Mie Theory for Scattering from Polydispersions 28 

3.3 Computational Methods 30 

3.4 Properties of Clay Samples i 32 

3.4.1 Physical Characteristics of Clay Samples 32 

3.4.2 Particle Size Distributions 35 

3.5 Results of Computations 45 

3.5.1 Volume Scattering Functions 47 

3.5.2 Volume Scattering Distribution Functions 54 

4.0 DEPENDENCE OF UPWELLING RADIANCE ON SCATTERING 

FUNCTIONS 61 

4.1 Results 63 

APPENDIX A RADIATIVE TRANSFER COMPUTER PROGRAM 71 

RADIATIVE TRANSFER COMPUTER PROGRAM - Code 1 73 

RADIATIVE TRANSFER COMPUTER PROGRAM - Code 2 89 

APPENDIX B MEASUREMENT OF SCATTERING FUNCTIONS 103 

B.l Scattering Functions 103 

B.2 Determination of Volume Scattering Function 103 

B.3 Scatterance Meters 104 

B.3.1 General Angle Scattering Meter 105 

B.3. 2 Small Angle Scattering Meter 106 

V 



APPENDIX C 

APPENDIX D 

APPENDIX E 

REFERENCES 

DISTRIBUTION 


LISTINGS FOR POLYMIE AND DBMIE ROUTINES 
USED TO CALCULATE THE VOLUME SCATTERING 
FUNCTIONS 

PROGRAM LISTING FOR CURFIT ROUTINE USED 
TO FIT THE THEORETICAL SIZE DISTRIBUTIONS 
TO THE EMPERICAL DATA 

RELATIONSHIP BETWEEN EXTINCTION, SCATTERING, 
AND ABSORPTION COEFFICIENTS AND THE MIE 
PARAMETERS 


LIST 


109 

123 

137 

141 

143 



LIST OF ILLUSTRATIONS 


Figure Number page 


1“1 Schematic Diagram o£ LaRC's Experimental 

Set-Up 3 

2-1 Scattering Functions Observed In-Situ 12 

2-2 Scattering Functions Observed in Deep 

Clear Oceanic Water and Very Turbrd 
Harbor Water 13 

2-3 Scatterxng Functions Observed In- Vitro 15 

2-4 Scattering Functions at Different 

Wavelengths 17 

2-5 Scattering Probability Functions Obtained 

in Deep Clear Oceanic Water and Very 
Turbid Harbor Waters 18 

2-6 The Effect of Scattering and Absorbing 

Materials in Fresh Water on Scattering 
Probability Function 19 

2-7 Scattering Probability Function Measured 

in Sargasso Sea at 2 Wavelengths 21 

2- 8 Scattering Probability Functions For 

Natural Ocean Waters and for Fresh 

Water that has Been Filtered and 

Artificially Modified 22 

3- 1 Cumulative Size Distribution of Feldspar 

Sample 36 

3-2 Cumulative Size Distribution of Calvert 

Sample 37 

3-3 Cumulative Size Distribution of Ball 

Sample 38 

3-4 Cumulative Size Distribution of Jordan 

Sample 39 

3-5 Particle Size Density Function for Feldspar 

(Modified Gamma Distribution) 41 

3-6 Cumulative Size Distribution Fit of Feldspar 

Sample Using Modified Gamma Distribution 42 

3-7 Particle Size Density Function for Ball Clay 

(Junges Distribution) 43 

3-8 Cumulative Size Distribution Fit of Ball Clay 

Sample Using Junge Distribution 44 

3-9 Volume Scattering Functions for Feldspar 

(A = 500NM) 48 

3-10 Volume Scattering Functions for Ball Clay 

(A = 500NM) 49 

3-11 Volume Scattering Functions for Feldspar 

(lOyM Cutoff A = 500NM) 50 

vii 



Figure Number 


Page 


3-12 

3-13 

3-14 

3-15 

3-16 

3-17 

3-18 

3-19 

3- 20 

4- 1 
4-2 
4-3 


Volume Scattering Functions for Ball Clay 
(lOyM Cutoff X = 500NM) 

Volume Scattering Functions for Feldspar 
(lOyM Cutoff X = 600NM) 

Volume Scattering Functions for Ball Clay 
(lOyM Cutoff X « 600NM) 

Volume Scattering Distribution Functions 
for Feldspar (X = 500NM) 

Volume Scattering Distribution Functions 
for Ball Clay (X = 500NM) 

Volume Scattering Distribution Functions 
for Feldspar (lOpM Cutoff X = 500NM) 
Volume Scattering Distribution Functions 
for Ball Clay (lOyM Cutoff X = 500NM) 
Wavelength Dependence of Volume Scattering 
Distribution Functions for Feldspar 
(lOyM. Cutoff) 

Wavelength Dependence of Volume Scattering 
Functions for Ball Clay (lOyM Cutoff) 
Backscattered Radiance vs. Upper Limit 
of the Exit Angle for s = 6 Meter^ 


Backscattered Radiance vs. Upper Limit 
of the Exit Angle for s = 12 Meter”^ 


Ratio of the Backscattered Radiance for 


Upper and Lower Bounding Scattering 
Functions 


51 

52 

53 

55 

56 

57 

58 

59 

60 
66 

'67 

69 


viii 



LIST OF TABLES 


Table Number 


3.1 Chemical Composition and Index of Refraction 

of Clay Constituents 

3.2 Settling Velocities of Sand and Silt In 

Still Water 

4.1 Optical Parameters Used in the Backs cat tered 

Radiance Calculations 


Page 

34 

46 

64 


IX 



1*0 INTRODUCTION AND CONCLUSIONS 


The importance of continuous monitoring of environmental water 
quality has long been recognized* The recent emphases placed on 
such operations are due to newly gained insights (1) in the 
limitations of the cleansing capability of the natural waters, (2) a 
better understanding of ecological consequences of water pollutants, 
and (3) availability of better information for assessing economic 
impacts of various stresses imposed on the water systems. Considering 
the dynamic character of the environmental waters the monitoring 
procedures for measuring water quality parameters should be based on 
timely data collection systems, such as can be provided by applica- 
tions of remote sensing technology. 

Hypothetically, in a remote sensing experiment the optical 
sensor measures the radiance signal which contains information on 
spectral and spatial variation of the source of radiation and the 
intervening media* The received radiance is then "processed" 
according to an established scheme which is a quantitative analytical 
procedure, and the radiance characteristics are ultimately related 
to the desired parameters. 

The data interpretation techniques for remote measurement of 
water quality parameters are presently in preliminary stages. 

Although some attempts have been made to develop analytical 
procedures for data processing, a generally accepted processing 
scheme has not emerged. 



Among the quantities that effect the radiance characteristics 
measured by a remote sensing instrument are: 

* Atmospheric path radiances and signal transmission effects 

* Spatial and spectral variability of atmospheric constitu- 
ents such as particulates and molecular species 

* Sun angle 

* Characteristics of air-water interface 

* Vertical non-homo genity of water bodies and bottom reflection 
properties 

Considering these effects and the fact that aquatic environments 
change continuously with the complex interactions between wind, water 
and land masses; the development of data interpretation schemes, in 
support of remote sensing, necessitates field experiments and 
controlled laboratory experiments as well as radiative transfer 
modeling approaches. A variety of field experiments from low and 
high flying aircraft and from satellite platforms have been con- 
ducted, or planned for the immediate future. 

A laboratory program is presently being pursued at the NASA- 
Langley Research Center (LaRC) . The purpose of this program is to 
investigate the remote sensing of water quality parameters under 
controlled conditions. During the first phase of this program, 
remote sensing applications of suspended particulates (various types 
of clays) have been investigated. A schematic diagram of LaRC*s 
experimental set-up is shown in Figure 1-1. In this experiment, the 
beam of a solar simulator is deflected to illuminate a large water 
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tank filled with turbid water; the water turbidity in the tank is 
caused by the stepwise introduction of specific amounts of particu- 
lates. An overhead detector system including a spectrometer, elec- 
tronics, and a camera, measures the strengths and the characteristics 
of the upwelling radiance signal. 

In order to analyze the experimental results and to optimize 

the experimental conditions a radiative transfer model has been 

developed for the LaRC’s experimental arrangement at The METREK 

Division of The MITRE Corporation. The description of the modeling 

approach and the results of a sensitivity study concerning the 

optimized spot size to be illuminated by the solar simulator have 

(1 2 ) 

been reported in two earlier documents. * 

The present report deals with variations in the characteristics 
of the backscattered radiance as a result of changes in the scatter- 
ing function, for various waters. The scattering function represents 
one of the important optical parameters of the turbid water and various 
scattering functions may represent various types of turbid waters. 

In Section 1.3 specific goals of the present report are described in 
more detail. Before this is dona however, it is necessary to 
summarize some background material on (1) optical parameters of turbid 
waters, and (2) on our radiative transfer modeling procedure. These 
background materials are treated in Sections 1.1, and 1.2 respectively. 
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1.1 optical Parameters of Turbid Water 


In the absence of polarization the following parameters are 
necessary for optical characterization of turbid water medium: 

• total absorption coefficient, a. 

• total scattering coefficient, s. 

These coefficients have the dimension of meter The attenuation 
coefficient, a, is the sum of absorption and scattering coefficients. 
The last parameter of interest is the scattering function, a(0). 

This function specifies the angular pattern of the scattering of a 
collimated beam from an infinitesimal volume of turbid water. The 



1.2 Radiative Transfer Model 


The development of METREK^s radiative transfer model is based 
on a two step process which is described in this section. The 
adopted modeling procedure is geared toward handling turbid type 
waters, and toward saving the computer time necessary for model 
execution* The model development includes the following steps: 

Step 1 * The outgoing radiance distribution just above 
the air-water interface, due to a narrow beam 
transmission in the turbid water media is 
calculated using Monte Carlo simulation tech- 
niques . 
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.Step 2 . The outgoing radiance emerging from the area 
.within the detector's field-of-view, and 
'traveling in a direction coincident with the 
range of the detector'''s acceptance angle is 
calculated using the interface radiance 
distribution (Step 1) and Integrating over 
'the incident ^eam area. 

The advantage of this approach as compared to conventional 
Monte Carlo simulation approaches is that the narrow beam considera- 
tion .allows the production of a better set of statistics .within 
reasonable computer resources. 

1 . 2.. 1 Monte Carlo Simulation for .Narrow Beam Transmission 

The advances in laser technology in the last decade have led to 
a 'variety of theoretical considerations nf the narrow beam trans- 
jmission in the water media. In general, the theoretical approaches 
may ibe -sub-divided in two categori'es,, ’(1) analytical solution of the 
.equa'tion of radiative transfer and (2) Monte Carlo 'simulation methods. 

The Monte 'Carlo .simulation methods avoid many of the mathematical 
complexities involved in the analytical solution approach, and for 
this reason are more .appropriate for calculating the narrow beam 
transmission. This is even more true in calculations simulating 
laboratory experiments where the experimental conditions, such as 
the tank ;geometry, significantly complicate the boundary conditions 
for the solution of the radiative .transfer equation. Thus, the 
Monte Carlo simulation method has been used in the development of 


'6 




the analytical model for LaRC’s experxment* A description of the 

Monte Carlo simulation approach, which is geared toward decreasing 

computer time and handling turbid rather than oceanic type waters 

is given in References 1 and 2. The procedure leading to the 

calculation of radiance is based on making use of the distribution 

of the emerging photons generated by the Monte Carlo program, 

(1 2 ) 

and the geometry of LaRC’s experimental arrangement* ^ 

The listing of the complete computer program, description of 
the input data, output data, and instructions for analysis of the 
output data to arrive at the upwelling radiance, are given in 
Appendix A. 

1, 3 Conclusions and Organization of This Report 

(1 2 ) 

In our previous reports, ^ we have documented the results 

of our modeling effort concerning the relationship between the spot 

size of the incident beam and the upwelling radiance, in the LaRC’s 

laboratory experiment* These results, however, were based on the 

(3) 

usage of only Morrison’s scattering probability function. In the 

present work we report on the effects of various inputs of both 
measured and calculated scattering probability functions. 

In Sections 2 and 3 we have (1) summarized the available 
information on the measurements of the scattering function, and (2) 
have utilized the Mie formalism to calculate the scattering function 
for polydispersed suspensions on the basis of size distribution 
measurements provided through the LaRC laboratory program, and 
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reasonable inputs for the index of refraction including its imaginary 

part. The compiled measured scattering probability functions for 

natural water, cover a wide range of turbid waters and show 

considerable variations. The upper and lower bounding measured for 

the scattering probability functions correspond to San Diego Harbor, 

and sea water filtered thoroughly. The scattering probability 

(3) 

function measured by Morrison^ , used in Reference (1,2) lies 
between these limits, closer to the upper bound. Due to the lack 
of sufficient observations no conclusions could be drawn regarding 
the changes of the measured scattering functions with wavelength. 

The calculated results of the scattering probability functions 
have been obtained for the following cases and their combinations: 

, Size distributions including large particle sizes (‘'^100 y) 

, Size distributions including a cutoff at 10 y 

• Zero or 0.004 for the imaginary part of the index of 
refraction 

• Two wavelengths values at 500 and 600 nm 

The conclusions derived from these results are : 

1) Size distributions including large particles sizes 
('vlOO y) lead to an extremely large forward scattering 
peak, which shows up as a fast rise in the scattering 
probability function. The scattering probability 
function calculated for this situation is higher than 
the upper bound of the measured functions as may be 
seen by comparing Figures 3-15 and 2-8. 

2) Size distributions including a cutoff at 10 y results 
in the scattering probability fmctions which lie 
between the upper and lower bounding of the measured 
probability functions shown in Figure 2-8. 
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3) The effect of non-zero imaginary part for the index of 
refraction is to decrease the fast rise of the scattering 
probability function at small angles, and to put these 
functions within the bounds of the measured data. 

4) The functions calculated for wavelengths of 500 nm and 
600 nm do not show significant differences. 

Based on the results and the conclusion described above three 
functions were selected for the investigation of the dependence of 
the upwelling radiance on the scattering function. These functions, 
which have been used in the Monte Carlo simulations radiative 
transfer code of Appendix A are: 

• The lower bound of the measured scattering probability 
function 

• The upper bound of the measured scattering probability 
function, and 

• The upper bound of the calculated scattering probability 
functions. This function has been calculated for Feldspar 
soil, a zero value for the imaginary part of the index of 
('ilOO p). This function is higher than the upper bound 

of the measured scattering functions. 


The turbidity levels treated in section 4.0 correspond to 
scattering coefficients s = 6.0 and s = 12 meter the wavelength 
of interest is 500 nm. The maximum number of incident photons 
traced in most computer runs is 10,000. The values calculated with 
the input of calculated upper bound scattering function is in good 
agreement with the measured upper bound scattering function for 
larger range of exit angle. However, for small range of exit angle 
(< 25^ degrees for s = 6.0 meter ^ and < 35^ for s = 12 meter no 
statistically significant result could be drived for this function, 
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from the ensemble of backscattered events for 10,000 incident photons. 
For thxs reason only the results derived from the use of measured 
upper and lower scattering functions were processed further, and 
form the basis of our conclusions. 

_X 

The results generated for s = 6.0 and s = 12.0 meter are in 
very good agreement as shown in Figure 4-3 (the figure-of-merit) , 
where the ratio of the backscattered radiances (radiance due to the 
lower limit measured scattering probability function, divided by the 
radiance due to the upper liniit measured scattering probability 
function) have been displayed as a function of the upper limit of the 
exit angle. As can be seen from Figure 4-3, the influence of the 
scattering probability function is quite significant, but decreases 
with decreasing exit angle. We expect that this trend will continue 
to be true for smaller angular ranges (such as 0 to 0.5° which 
represent the acceptance angle of the LaRC's overhead detector) and, 
therefore, conclude that the effect of various scattering probability 
functions is not significant in the LaRC's experimental set-up. 

The reason the smaller angular ranges were not examined specif- 
ically in this report has been due to the constraint on computer 
resources. It is recommended, therefore, that the computer program 
developed in this report be executed for a larger number of photons 
(larger than 10,000 photons considered in this study) to strengthen 
our conclusions. 
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2.0 MEASURED SCATTERING FUNCTIONS 


In. order to measure the complete scattering function, the scat- 
terance must be observed at a number of angles between 0° and 180°. 
Two types of scattering meters have been used in the past for the 
measurements of the scattering functions. These are: (1) general 

angle scatterxng meters, and (2) low angle scattering meters. The 
mathematical definition of scattering function and an overview of 
the scattering meters are given in Appendix B. 

The instrumentation required for in-situ measurements of the 
scattering functions are very sophisticated, hence only a small 
number of such measurements have been performed. 

Figures 2-1 and 2-2 show several in-situ measured scattering 
functions covering turbid to clear water conditions. Figure 2-1 
represents the observations made in lake water, coastal waters, 
Atlantic surface water, Pacific near-coastal water, Mediter- 
ranean, and Saragasso Sea water. Most of these observations 
are taken between 0 = 10° and 6 = 155°. Figure 2-2 illustrates the 
measurements taken by the Scripps Institution of Oceanography in 
deep clear oceanic water (tongue of the ocean), near shore ocean 
water (off shore of Southern California), and very turbid harbor 
water (San Diego Harbor). The measurements shown in Figure 2-2 are 
\ 'cArffed/oub-iO.ver the entire range 0°< 9 < 180°. The scattering 
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FIGURE 2-2 

SCATTERING FUNCTIONS OBSERVED IN DEEP CLEAR OCEANIC 
WATER AND VERY TURBID HARBOR WATER 





functions shown in both these figures are similar in general form* 

The differences between the scattering functions are most pronounced 
in the backward region above 90^, and in the forward region below 
IQo. Although the differences do reflect real variatxons in the 
scattering functions for various areas, they may reflect the inherent 
experimental difficulties as well* 

The experimental difficulties become more striking when the 
scattering functions are measured in-vitro* The observations taken 
in-vitro by Petzold^^^^ are probably the most reliable ones* The 
measurements were taken to determine the effect of adding scattering 
and absorbing materials in the water* For this, scattering materials 
(compounds of aluminum hydroxide and magnesium hydroxide), and 
absorbing materials (black dye nigrosin), were introduced into fresh 
water pumped through a filter containing diatomaceous earth* The 
resultant change in scattering functions as observed with the 
scattering meters are presented in Figure 2-3* It is clear from 
Figure 2-3 that the scattering functions are insensitive to the 
absorption properties of the water. 

2.1 Variation of Scattering Function with Wavelength A 

Not many of the experiments either in-situ or in-vitro so far 
have been performed for different wavelengths. Most of the obser- 
vations are in 460 ^ \ 4 655 nm wavelength region. The scattering 
functions presented in Figures 2-2 and 2-3 were measured at 
A = 530 nm. Due to a lack of observations at other wavelengths for 
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the same meters and under similar conditions (Figures 2-2 and 2-3) , 
it is difficult to draw any conclusions regarding the changes in 
scattering functions with wavelength ♦ However, Kullenberg has 
measured a(0) at 655 nm and 460 nm in the Sargasso Sea* These 

measurements are shown in Figure 2-4. The scattering function is 
evidently the same at both these wavelengths, in the forward scat- 
tering region of 0 ^ 0 ^ 35^. 

2.2 Scattering Probability Function F(6) 

The scattering probability fimction, F(9), has been defined by 
equation (1-1). F(0) is the ratio of power scattered into angles less 

than 0 relative to the total power scattered in all directions. F(0) 
is an important parameter and is a measure of forward as well as back- 
ward scattering in natural environment waters. F(0) is the function 
used in the Monte Carlo simulation model, as mentioned in the intro- 
duction. 

Figure 2-5 shows the scattering probability function obtained 
by integrating the function represented in Figure 2-2 , while Figure 
2-6 illustrates F(0) obtained from Figure 2-3. The probability 
scattering functions presented in Figure 2-6 show the effect of 
adding scattering and absorbing materials in the waters. Clearly, 
the addition of scattering material increases the backs cat tering 
whereas addition of absorbing material contributes insignificant 
changes to the scattering probability function. 
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Scattering Function o(0), la str 


Scattering Function Measured in 
Sargasso Sea at 460 nin. 

Scattering Function Measured in 
Sargasso Sea at 655 mn. 
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FIGURE 2-4 

SCATTERING FUNCTIONS AT DIFFERENT WAVELENGTHS 
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FIGURE 2-5 

SCATTERING PROBABILITY FUNCTIONS OBTAINED IN DEEP CLEAR OCEANIC 
WATER AND VERY TURBID HARBOR WATERS 




FIGURE 2-6 

THE EFFECT OF SCATTERING AND ABSORBING MATERIALS IN FRESH WATER 
ON SCATTERING PROBABILITY FUNCTION 



To illustrate the effect of varying wavelength on the scattering 
probability fvinction, the functions presented in Figure 2-4 were 
filled in the angular range larger than 1° degree, were extrapolated 
into the angular range smaller than 1° degree, and were integrated 
to obtaxn F(9) at 655 nm and 540 nm, as shown in Figure 2-7. 

Since one of the objectives of this paper is to record experi- 
mentally determined upper and lower bounding scattering probability 
functions, the information from Figures 2—5, 2—6, and 2—7 are shown 
collectively in Figure 2-8, The lowest bound on the scattering 
probability function is given by the pure water, where particulate 
scattering is insignificant. Natural environmental waters are not 
usually free of particulates and, therefore, experiments have been 
performed to define their characteristics. An experiment conducted 
at Scripps Institution of Oceanography examined sea water pumped 
into the laboratory and measured scattering probability functions 
for the water as delivered, and after several steps of filtration. 
After 18 hours of filtering, low-angle foward scattering signals 
have been found too low to be measurable. The results obtained from 

this experiment, in addition to the scattering probability ftinctions 

(3) 

obtained by Morrison at Long Island Sound stations, are included 
in Figure 2-8. 

The work presented in this section indicates that the San Diego 
Harbor water, the most turbid water, gives the upper bound to the 
experimentally determined scattering probability functions. The 


20 



ORIGINAL PAGE IS 
OP POOR QUAIZIY 



FIGURE 2-7 

SCATTERIMG PROBABILITY FUNCTION MEASURED IN SARGASSO SEA 

AT 2 WAVELENGTHS 




FIGURE 2-8 

SCATTERING PROBABILITY FUNCTIONS FOR NATURAL OCEAN WATERS 
AND FOR FRESH WATER THAT HAS BEEN FILTERED 
AND ARTIFICIALLY MODIFIED 




lower bound is given by the sea water thoroughly filtered. The 
scatterance characteristics of various waters considered are quxte 
different. Very turbid waters show very high forward scatterance. 

At 1° scattering angle, the forward scattering measured in San Diego 
Harbor water is almost 15 times of that measured in filtered water. 
This ratio reduces to three at 10° scattering angle. 

The implications of these results on remote detection of water 
turbidity will be discussed in Section 4.0. 
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3*0 CALCULATED SCATTERING FUNCTIONS 

This section is devoted to the theoretical treatment of 
scattering and absorption from suspended particulates. The Mie 
theory of light scattering from a single particule is treated in 
Sub-section 3.1. The extension of Mie theory to the case of 
polydispersed suspensions is then discussed along with the compu- 
tational methods used to calculate the scattering function » in Sub- 
sections 3.2 and 3.3 (Appendix E discusses the relationships between 
the Mie parameters and the extiction, scattering and absorption 
coefficients). Sub-section 3,4 includes a discussion of the size 
distributions and optical properties of the clay sediments 
considered in the calculations. Finally, in Section 3.5 the results 
of the calculation of the scattering function are presented along 
with a discussion of their implications for the NASA/Langley tank 
experiment. 

The following discussion of the Mie theory of scattering and 
the computational methods is a brief summary. For more detailed 
discussions of Mie theory for single scattering, the reader is 
referred to References 13, 14, and 15. Reference 16 contains a 
good discussion of Mie scattering from polydispersions and reference 
17 contains the details of the computational procedures and 
requirements • 
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3.1 Mie Theory for Single Particle Scattering 


When light is incident on a particle, it undergoes both scat- 
terxng and absorption (we will ignore inelastic scattering processes 
which result in a change in frequency). The characteristics of the 
scattered radiation depend on the wavelength, X, of the incident 
light, the generally complex index of refraction, m, of the particle, 
and the size, r, and shape of the particle. In this report we will 
restrict the discussion to spherical particles; for the treatment of 
inorganic sediments in water this is probably not a serious 
restrictxon. 

If a monochromatic beam of light of xntensity is xncident on 
a spherical particle at an angle 0=0, then the scattered intensity 
xs given by 


I(x,m, 0) 



a(x,m,0) 


(3-1) 


tfliere a (x,m, 0) is the single particle scattering function, a (x,m, 6) 
depends in general, on the size parameter. 


X 


2Trr 

X 


(3-2) 


and the complex index of refraction, m. The calculation of a (x,m, 6)’ 
requxres the solution of Maxwell’s equation in spherxcal coordxnates 
with a discontinuous change in the index of refraction across the 
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spherical surface. This solution was originally derived by G. 

(19') 

and independently by P. Debye^ \ 

The scattering function can be written as: 


a(x,m, 0) 



(x,m,e) + 
_ 


(x,m, e) 


and the Mie solution is 


a^(x,m, e) 
a^Cx.m, e) 


S^(x,m,0) (x,m, 0) 

* 

S2(x,m,0) (x,m,e) 


(3-3) 


(3-4) 


Where S^(x,m, 9) and S2(x,m, 6) are the complex amplitudes for the 
scattered radiation, 


S^(x,m,8) 


“ (2n+l) 

n(n+l) 


S„(x,m, 0) = Z 

n=l 


(2n.+l) 

n(n+l) 


a^(x,m)Tr^(p)+b^(x,m)x^(y) 


I 


L 

|b^(x,m)ir^(u)+a^(x,m)T^(y) „ 


(3-5) 


In these expressions derivatives of the Legendre 

polynomials : 






dP^(y) 

dy 


yiT^(y) - (1-y ) 


dir (y) 
n 

dy 


(3-6) 
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(where \x = cos 0 ). Also 


a^(x,m) if) ' (rax) E (x) - nn|> (mx) E " (x) 

n n n n ^^- 7 ) 

ini|)^(inx)ii)^(x) - ij>^(rax)<|)'(x) 
b^(x,m) - m^^(njx)5^(x) - ij;^(nix)E^(x) 


and the tj;'s and E's are related to the spherical Bessel functions of 
the first and second kinds (j^ and respectively) : 


En<x) = X j^(x)-iy^(x) 

5'(x) = 


3.2 Mie Theory for Scattering from Poly dispersions 

A polydispersion is a suspension of scattering particles of 
uniform physical characteristics but of varying number concentration 
depending on particle size. Because of the existence of different 
particle sizes it makes little sense to talk of scattering from a 
single particle. Instead, it is useful to consider the scattering 
properties of a small volume element containing a number of particles. 
The size of this volume element is of some, at least theoretical, 
importance. Clearly, If it is to be used to represent the scattering 
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properties of all similar volume elements then it must contain a 
representative set of particle sizes - this requires that the volume 
element not be too small* On the other hand, since we are consider- 
ing only single scattering from the volume element, it must not be 
too large. An additional condition that must be imposed is that the 
inter-particle separation be large compared to the wavelength. The 
reason for this is that the interaction of light with a particle 
will be assumed independent of the interactions with all other 
particles. This condition requires that the particle density in the 
volume element not be too large. For our purposes, it will be assumed 
that all of the above conditions are satisfied. 

The polydispersion can be completely specified, for our purposes, 
by an index of refraction m and a probability density function n(r). 
The density function gives the relative concentration of each size 
contained in a volume element. 

The characteristics of the scattered radiation due to the volume 
element can then be represented by a volume scattering function 
a(m, 0) in a manner analogous to Equation (3-1): 

I(m,e) = — ^ a(m,0)I^ (3-9) 

Air 

The scattering function can be calculated from the set of partxcle 
scattering functions: 


a(m, 6) 



( 3 - 10 ) 
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where 


J n(r)dr = N (3-11) 

o 

and N is the total number of particles per unit volume. In what 
follows, N will be assumed to be unity since a(m, 6) scales with N, 

The ability to represent a(m, 0) as a linear superposition of the 
a(s,m, 0)s is a direct consequence of our assumption that the inter- 
particle separation is much greater than 

The calculation of a(m,9) thus reduces to calculations of the 
individual a(s,m,0) and then integration over all sizes with the 
proper weighting given by n(r). 

3.3 Computational Methods 

The calculation of the scattering functions and the averaging 
over size distributions was carried out on an IBM 370/148. The 
program listings are reproduced in Appendix C. 

In computing the sums in Equation (3-5), the major difficulty 
arises in the evaluation of the a^(x,ra) and b^(x,m). Using the 
definitions of i{)^, 5^, and and the standard recurrence 

relations for the -Bessel functions. Equation (3-7) can be rewritten: 



rA (mx) ■] 

|- Re 


-Re 


W <mx) 1 




|^inA^(mx)+n/xj 

1- Re[5„(x)] 

-Re[c^-l(x) 


finA (rax)-Hi/x^ 

L " 







Where 


A (mx) 
n 


ip^(mx) 


(3-13) 


the logarithmic derivative of iJ<^(Tnx) , and Re denotes the real part. 

The natural approach to the evaluation of Equation (3-12) is to employ 
a standard upward recurrence procedure. Unfortunately, if the 
imaginary part of m, Ira(ra), is not zero and n is larger then the up- 
ward recurrence procedure results in larger instabilities in the 
calculation of A^(mx). For this reason, the DBMIE subroutine employs 
a downward recurrence procedure to calculate the A^(rax)s. These 
values are then stored for use in the evaluation of Equation (3-12) . 
Because of the large storage requirements resulting from this proce- 
dure (n 7000) , and the fact that double precision is employed in 
all of the calculations, a virtual machine with 512 K bytes of 
storage is required for the implementation of the DBMIE and POLYMIE 
routines . 


While the scattering functions are computed in the DBMIE 
subroutine, the average, Equation (3-10), is computed in the calling 
routine POLYMIE. While analytic functions have been used for the 
Size distributions, n(r), the integral has been approximated by a 
summation over a discrete set of radii. Tests to determine the 


effect of using a summing procedure have shown that this results in 
no loss of accuracy. In addition, test runs were made to compare 
the results when Ar = O.ly (0,1 micron) and Ar = ly -were used in the 
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summing procedure. The use of Ar =? ly resulted in no significant 
change in the results from those obtained using Ar = O.lp over the 
range 0 <-r < lOOy. Calculations were made using r = lOOp 

TQ3-X 

(Ar ^ Ip) and r - lOp (Ar ?= O.lp). A discussion of the proper upper 
max 

limit for r is given in Section 3.4* 

The amount of virtual CPU time required for these calculations 
is significant and has been a major factor in determining 3.nd 

Ar. As an example > the calculation of the volume scattering function 
for a polydispersion with m = 1.144 - O.Oi, X =0.5]i, r - lOOy 

and Ar = Ip requires approximately 26 minutes of virtual CPU 
time. 

3.4 Properties of Clay Samples 

Data on four different clay samples were provided by NASA/LaRC. 
This data consisted of empirical size distribution curves as well as 
brief descriptions of chemical composition. The physical character- 
istics of the clay are discussed in Section 3.4.1 while the size dis- 
tributions are presented in Section 3.4.2. 

3.4.1 Physical Characteristics of Clay Samples 

Four types of clay were selected by NASA/LaRC. These were: 
Feldspar, Calvert, Ball and Jordan. According to the analysis of 
these clays performed by NASA/LaRC the compositions are: 

* Feldspar - Feldspar and Quartz minerals 

• Calvert and Jordan - Kaolinite and Illite 

, Ball - Montmorilloite, Kaolinite and Illite 
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The real refractive index and chemical components of these 

f 21') 

minerals is shown in Table 3.1 . For reasons which will be dis- 

cussed in Section 3.4.2, Feldspar and Ball clay were chosen to be 
included in this study. 

To estimate the index of refraction of the clay samples, we 

take a simple average of the indices of refraction of the components. 

Thus, for both Feldspar and Ball clay, the real part of the index of 

refraction is estimated as 

RE (ra,. ) = 1.53 
Axr^ 

This, of course, is the index of refraction with respect to air and 
we require the index of refraction with respect to water which can 
be obtained by dividing t>y the index of refraction of water 

1.337 (for wavelengths of approximately 500 nra) . 

Thus 


Re(m „ ) = 1.144 

water 


Estimating the imaginary part of the index of refraction is not 

so straightforward, since direct measurements of Im(m) have not been 

made. Since these minerals have very low conductivity, it is expected 

that the imaginary part of m will be quite small. The imaginary part 

of m has been measured for soil aerosols and has been found to be 

( 22 ) 


about .005 (with respect to air), 
Im(m) will be used: 


For- this study two values for 


Im(m ^ ) = 

water 


0 , Non- absorbing 

0.005 


L1.337 


= 0.004, Weakly-absorbing 
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TABLE 3.1 


CHEMICAL COlffOSITION AND INDEX OF EEFRACTION 
OF CLAY CONSTITUENTS 


NAME 


Kaolinite 


CHEMICAL COMPOSITION 
Al2O3.2SiO2.2H2O 


INDEX OF REFRACTION 
1.56 


Illite 3 Al 2 ^Sx^_g 3 A 1 ^_^^ 302 ( 0 H)^ 1.54 

Montmorilloite (.5Ca,Na) ^)Al,Mh,Fe)^(sx,Al)g)2Q(H0)^nH20 1.48 


Feldspars'; 

Microcline 


Andesine 


Anthoclase 


K20.A1203.6S102 
(Ca 03 ^Na 20 ) AI 2 O 3 . 4S102 
(Na,K)20,Al203.6Si02 


1.52 
1.55 

1.53 
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3.4t2 Particle Size Distributions 

Emperical cumulative size distributions for the four sairples 
were provided by NASA/LaRC and are sho™ in Figures 3-1 > 3-2, 3-3 and 
3-4. It is apparent from these figures that the size distributions 
for Ball, Jordan, and Calvert differ significantly from the size 
distribution for Feldspar. Since it was planned that two distri- 
butions would be employed, Feldspar and Ball clay were chosen. 

This choice allows the investigation of the effect of radically 
different size distributions. 

To utilize the size distribution information, it is necessary 

to determine size distribution density functions, n(r), which specify 

the relative number of particles with radius r per unit volume. If 

we denote the cumulative size distribution as provided by NASA/LaRC 

as N(r ) then the relationship between N(r ) and n(r) is given by: 
o o 


or 


N(r^) = 1- n(r) dr, 

0 


, ^ dN(r ) 
n(r) - -jj o 

o 



(3-14) 


(3-15) 


A general curve fitting routine (See Appendix D) was used to deter- 
mine the best distribution for both the Ball clay and Feldspar. 

For the Feldspar sample, it was found that the data was well 
represented by a modified Gamma distribution: 

, , ^2 I (3-16) 

n(r) = a^r exp l-a^r 
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Percentage of Particles 



Particle Size, Milluneters 


FIGURE 3-1 

CUMULATIVE SIZE DISTRIBUTION OF FELDSPAR SAMPLE 



Percentage ox 



FIGURE 3-2 

CUMULATIVE SIZE DISTRIBUTION OF CALVERT SAMPLE 





FIGURE 3-3 

CUMULATIVE SIZE DISTRIBUTION OF BALL SAMPLE 



Percentage 



Particle Sxze, Mxllimeters 


FIGURE 3-4 

CUMULATIVE SIZE DISTRIBUTION OF JORDAN SAMPLE 




The parameters were determined^ using a minimum mean square error 
criterion, to be 

- 2.05089 
^2 ' 0.671066 
= 3.58393 
= 0.218499 

A plot of this size distribution density function is shown in Figure 
3-5, while a plot of the corresponding cumulative size distribution 
function (as obtained from Equation 3-16) is shoxm in Figure 3-6. As 
can be seen in Figure 3-6, the modified Gamma distribution gives a 
good fit to the data points obtained in the NAS A/LaRC analysis. 

To fit the size distribution of the Ball clay sample, Junge^s 
distribution model was chosen; 

n(r) = r“^2 (3-17) 

With the parameters, 
a^ = .2006 

a^ = 1.624746 

determined using the same curve fitting routine employed for Feldspar. 
The size distribution density function and the cumulative size dis- 
tribution function for Ball clay using Junge’s distribution are shown 
in Figures 3-7 and 3-8. It is apparent from Figure 3-7 that Junge^s 
distribution function is not, strictly speaking, a probability dis- 
tribution since the integral (Equation 3-11), 

00 

J n(r) dr = N 
0 
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n.(r) ~ Relative Humber of Particles per Unit Volume 



Particle Size, Microns 


FIGURE 3-5 

PARTICLE SIZE DENSITY FUNCTION FOR FELDSPAR {MODIFIED GAMMA DISTRIBUTION) 
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Percentage of Particles 
Coarser by Weight 



.1 .05 .02 .01 .005 .002 .001 .0005 .0002 .0001 

Particle Size, Millimeters 


FIGURE 3-6 

CUMULATIVE SIZE DISTRIBUTION FIT OF FELDSPAR SAMPLE USING MODIFIED GAMMA DISTRIBUTION 



n(r) - Relative Number of Particles per Unit Volume 



FIGURE 3-7 

PARTICLE SIZE DENSITY FUNCTION FOR BALL CLAY {JUNGES DISTRIBUTION) 
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FIGURE 3-8 

T OF BALL CLAY SAMPLE USING JUNGE DISTRIBUT 


CUMULATIVE SIZE DISTRIBUTION FI 


ION 





can not be normalized, i.e., N is infinite. However, Junge's distri- 
bution has been found to accurately represent particle sizes of ocean 
(23) 

sediments. In addition, the lower and upper limits of integration 

in Equations (3-11) and (3-10) are not set equal to zero and infinity, 

in practice, allowing Equation (3-11) to be normalized. 

The question of the proper upper limit for Equation (3-10) and 

Equation (3-11) is of more than theoretical interest. From the 

empirical size distributions provided by NASA/LaRC, it appears that 

an upper limit in Equation (3-10) should be chosen as 100 microns 

(24) 

(ym). However, as can be seen in Table 3.2 the settling rate 
for 100 pm particles is on the order of forty seconds. Thus, the 
history of the particulates in the body of water is important. If 
the particulates have been allowed to settle, then the size distri- 
butions determined before the particles are introduced into the water 
are inappropriate. In the NASA/LaRC water tank experiment the water 
is continuously mixed, thus forcing the large particles to remain in 
suspension. In order to investigate the effect of settling, two 
upper limits, 100 pm and 10 pm, were chosen for the integrals of 
Equations (3-10) and (3-11) . Equation (3-11) was used to properly 
normalize Equation (3-10) with respect to the choice of upper limit, 
3.5 Results of Computations 

. . qV'*' ’* 

. ‘ 1 • ' ' ' ' 

The' fiesu'lts of the ’computation of the volume scattering 

4.1. , i ^ ' 

functions (3.5.1) and the volume scattering distribution functions 
(3.5.2), using the size distributions of Section 3.4, are presented 
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TABLE 3.2 


settling velocities 


OF SAND AND 


SILT IN STILL WATER 


(Source- Amer Water Works Assoc.) 

(Temperature SO°F, all particles assumed to have a specific gravity of 2 65) 


Diameter of 


Settling 

Time Required to 

partigte 

Order of Size 

Velocity 

Settle 1 Foot 


mm. 


mm /sec 


100 

Gravel 

1,000 

0 3 seconds 

IsO 


100 

3 0 seconds 

08 


83 


06 


63 


0 5 
04 

^ Coarse Sand 

53 

42 


03 ! 


32 


02 j 


21 


0.15 < 


15 


0 10 ' 


8 

38 0 seconds 

0 08 


6 


0 06 


38 


0 05 

\ Fine Sand 

2 9 


0 04 

2 1 


0 03 


1 3 


0 02 


0 62 


0 015 J 


0 35 


0 010 > 


0 154 

33 0 minutes 

0 008 


0 098 


0 006 


0 065 


0 005 

S Silt 

0 0385 


0 004 

0 0247 


0 003 


0 0138 


0 002 


0 0062 


0 0015 


0 0035 


0 001 

Bacteria 

0 00154 

55 0 hours 

0 0001 

Clay Particles 

0 0000154 

230 0 days 

0 00001 

Col'oidal Particles 

0 000000154 

63 0 years 
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in this section. In addition to examining the effect of settling 
on the calculations, the wavelength dependence of the scattering 
functions are also investigated. 

3.5.1 Volume Scattering Functions 

The computed volume scattering functions are shown in 
Figures 3.9 through 3.14. 

Figures 3.9 and 3.14 display the extremely large forward scat- 
tering peak which is primarily the result of including the large 
(o-lOO ym) particulates in the size distributions. Both the Feldspar 
and Ball clay phase functions show considerable difference between 
the non- absorbing and absorbing cases at large angle. While it is 
not evident in the figures, the forward scattering peak is larger 
for the absorbing case at small but non-zero angles (6~0.5°). 

Figures 3-11 and 3-12 demonstrate the effect of cutting the size 
distributions off at 10 ym instead of 100 ym. The relative size of 
the forward peak is reduced and the difference between the absorbing 
and non-absorbing cases at large angles is reduced. It is inter- 
esting to note that, although the shape of the Feldspar and Ball 
clay size distributions are very different, the upper limit on the 
size appears to be much more important in terms of the difference in 
phase functions. 

Figures 3-13 and 3-14 show the scattering functions computed for 

X = 600 nm (with a 10 ym cutoff) instead of X = 500 nm as in Figures 

3-11 and 3-12. It can be seen that the phase functions are not heavily 
'■ 1 ■ • u.dSi > 

’ ' ‘ I* ! * 1 ' J 
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Scattering Function 
(Un-normalized) 



FIGURE 3-10 

VOLUME SCATTERING FUNCTIONS FOR BALL CLAY (X = 500NM) 




Scattering Function 
(Un-normalized) 



FIGURE 3-11 

VOLUME SCATTERING FUNCTIONS FOR FELDSPAR (lO^tM CUTOFF X = 500NM) 
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FIGURE 3-12 

VOLUME SCATTERING FUNCTIONS FOR BALL CLAY (10/iM CUTOFF X = 500NM) 
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FIGURE 3-14 

VOLUME SCATTERING FUNCTIONS FOR BALL CLAY (lOjuM CUTOFF X = 600NM) 






wavelength dependent. In fact, it can be shown that for a uniform 
size distribution and upper and lower limits of zero and infinity in 
Equation (3--10) , the volume scattering function will be strictly 
independent of wavelength, 

3,5,2 Volume Scattering Distribution Functions 

While the volume scattering function describes the angular 
dependence of scattered radiation, a more important function for use 
with the Monte Carlo simulation is the volume scattering distri- 
bution function, F(0), defined by equation (1“1) , The distribution 
function gives the normalized probability that a photon is scattered 
in the range 0 to 0 degrees. The volume scattering distribution 
functions for the cases considered in Section 3,5 are shown in Figures 
3-15 through 3-20, 

It is again apparent in Figures 3-15 and 3-16 that there is a 
considerable difference between the absorbing and non-absorbing case. 

The difference due to the Feldspar and Ball clay size distributions 
is small. 

As with the scattering functions, the use of a 10 ym cutoff decreases 
the difference between the absorbing and non-absorbing cases. In 
addition, the volume scattering distribution functions are changed 
consdierably when the 10 ym cutoff is imposed. 

Figures 3-19 and 3-20 demonstrate the small change in the 
volume scattering distribution functions when the wavelength is 
changed. 
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FIGURE 3-15 

VOLUME SCATTERING DISTRIBUTION FUNCTIONS FOR FELDSPAR (X-500NM) 
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FIGURE 3-18 

VOLUME SCATTERING DISTRIBUTION FUNCTIONS FOR BALL CLAY (lO/xM CUTOFF X = 500NM) 
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FIGURE 3-20 

WAVELENGTH DEPENDENCE OF VOLUME SCATTERING FUNCTIONS 
FOR BALL CLAY CUTOFF) 




4,0 DEPENDENCE OF UPWELLING RADIANCE ON SCATTERING FUNCTION 

In this section we describe our results on the dependence of the 
upwelling radiance as it relates to the variations of the scattering 
function, or equivalently to its integrated form the scattering pro- 
bability function. Before this is done, however, we will summarize 
the information on the scattering probability functions derived earlier. 

In the previous two sections, we have (1) summarized the available 
information on the measurements of the scattering function, and (2) 
have utilized the Mie formalism to calculate the scattering function 
for polydispersed suspensions on the basis of size distribution 
measurements provided through the LaRC laboratory program. The 
compiled measured scattering probability functions for natural water, 
Figure 2-8, cover a wide range of turbid waters and show considerable 
variations. The upper and lower bounding measured for the scattering 
probability functions correspond to San Diego Harbor, sea water 
filtered thoroughly. The scattering probability function measured 
by Morrison \ used in Reference (1,2) lies between these limits, 
closer to the upper bound. Due to the lack of sufficient observations 
no conclusions could be drawn regarding the changes of the measured 
scattering functions with wavelenth. The calculated results of the 
scattering probability functions have been obtained for the following 
cases and their combinations: 

• Size distributions including large particle sizes (--'lOO \xm) 
m Size distributions including a cutoff at 10 pm 



• Zero or 0.004 for the imaginary part of the index of refraction 

• Two wavelengths values at 500 and 600 nm 

The conclusions which may be derived from these results are: 

1) Size distributions including large particles sizes 
(^100 u) lead to an extremely large forward scattering 
peak, which shows up as a fast rise in the scattering 
probability function. The scattering probability 
function calculated for this situation is higher than 
the upper bound of the measured functions as may be 
seen by comparing Figures 3-15 and 2-8, 

2) Size distributions including a cutoff at 10 p results 
in the scattering probability functions which lie 
between the upper and lower bounding of the measured 
probability functions shown in Figure 2-8. 

3) The effect of non-zero imaginary part for the index of 
refraction is to decrease the fast rise of the scattering 
probability function at small angles, and to put these 
functions within the bounds of the measured data. 

4) The functions calculated for wavelengths of 500 nm and 
600 nm do not show significant differences. 

Based on the results and the conclusion described above three 
functions were selected for the investigation of the dependence of 
the upwelling radiance on the scattering function. These functions, 
which were input to the Monte Carlo simulations radiative transfer 
code of Appendix A, have been designated by SCATR 1, SCATR 2, and 
SCATR 3. SCATR 2 is the lower bound of the measured scattering 
probability function shown in Figure 2-8. SCATR 1 is the upper bound 
of the measured scattering probability function shown in Figure 2-8. 
SCATR 3 is the upper bound of the calculated scattering probability 
functions, and is shown in Figure 3-15. This function has been 
calculated for Felspar soil, a zero value for the imaginary part of 



xndex of refraction, a size distribution including large particles 
(^100 y) at 500 nm wavelength. 

4*1 Results 

Besides the parameters characterizing the cross sectional radius 
(1.2 meters) and the height (2.6 meters) of the LaRC cylindrical 
water tank, and the reflectivity of the tank walls (3.0 percent) the 
following input parameters are required for the model: 

(1) Total scattering coefficient s , 

(2) Total absorption coefficient a, 

(3) Scattering probability function. 

A fourth model input concerns the maximum number of photons to be 
traced in each computer run. 

The results presented in the remainder of this section refer to 
two turbidity levels which have been simulated in the model. These 
turbidity levels correspond to s = 6 meter and s = 12 meter ^ 
respectively. The wavelength considered is 500 nm. From the 
functional relationship between a/s ratio and the wavelength, reported 
in Reference 1, the value of a/s for particles at 500 nm is 0,27. 

Based on this value, absorption coefficients of 1.6 and 3*2 meter ^ 

have been calculated for s - 6 and s - 12 meter ^ respectively, and 

are shown in Table 4*1. 

On making use of the computer code documented in Appendix A the 
radiances emerging from the area within the field of view of the over- 
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TABLE 4.1 


OPTICAL PARAMETERS USED IN THE BACKSCATTERED 
RADIANCE CALCULATIONS 



TOTAL 

TOTAL 

TOTAL 

WAVELENGTH 

SCATTERING 

ABSORPTION 

ATTENUATION 

(nm) 

COEFFICIENT 

COEFFICIENT 

COEFFICIENT 


(meter“l) 

(meter-1) 

(meter"l) 


s 

a 

a 

500 

6 0 

1.6 

7.6 


12.0 

3-2 

15.2 


64 



head detector in the LaRC^s experimental arrangement,* and into the 
exit angles in the range 0^10 degrees, 0-20 degrees, 0-30 degrees, 
have been calculated. The results of these calculations in terms of 
the backseat tered radiance vs, the upper limit of the exit angle is 
shown in Figures 4-1 and 4-2 for s = 6.0 and s = 12.0 meter , Three 
scattering probability functions, namely the measured upper and 
lower bounding functions have been used. The model has been executed 
for 10,000 photons in each case. The values calculated with the 
input of calculated upper bounding scattering function is in good 
agreement (the shape of the respective curves) with the measured 
upper bounding scattering function for the large range of the exit 
angles. For the small range of the exit angles, (< 25^ degrees for 
s - 6 meter ^ and < 35^ for s - 12 meter no statistically 
significant result could be derived from the ensemble of backseat tered 
photons for 10,000 incident photons. For this reason the reminder 
of this report will discuss the results concerning the upper two 
curves in Figure 4-1 and 4-2. 

The presented results indicate that the upwelling radiance has 
a strong dependence on the scattering function used. This dependence 
seems to get less important with decreasing range of the exit angle. 

If the same trend continues to be true for smaller than 10^ angles 

* 

An area 2,5 cm in radius in the middle of the incident spot which 
IS about 30 cm in diameter (see Fxgure 1-1 for reference) . The 
incident beam impinges upon the water surface at an angle of 13.5 
degrees xn the air (9,0 degrees in the water)’. 

65 


original pag^ 
0 ^ POOE 



Backscattered Radiance 
(Arb. Units) 


1 


.1 



S = 6 meter 

Incident Spot Radius - 15 cm 
X = 500 nm 


J L 

10 20 

Exit Angle, Degrees 


FIGURE 4-1 

BACKSCATTERED RADIANCE VS. UPPER I 
THE EXIT ANGLE FOR s = 6 METER 
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FIGURE 4-2 

BACKSCATTERED RADIANCE VS. UPPER LIMIT OF 
THE EXIT ANGLE FOR s = 12 METER-”* 
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than 10° angles (for which no significant statistic could be derived 
for 10,000 photons)* then, at 0.5° angle which is the actual accep- 
tance angle of the LaRG's detector, the effect of various scattering 
functions will not be significant. This is displayed graphically 
by the results presented in Figure 4-3, where the ratio of the 
backscattered radiances for the upper and lower bounding of the 
scattering function is shown as a function of the upper limit of 
the exit angle. These results will be discussed in more detail in 
section 1.4. 


2D, 000 photons were traced to produce the results shown in the 
lower curve in Figure 4-2. 
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FIGURE 4-3 

RATIO OF THE BACKSCATTERED RADIANCE FOR UPPER AND LOWER 
BOUNDING SCATTERING FUNCTIONS 
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APPENDIX A 


RADIATIVE TRANSFER COMPUTER PROGRAM 

In this appendix we have included two versions of our radiative 
transfer code. These programs are appropriately modified versions of 
the program listed in Reference 2. The modified conputer codes make 
it easier to incorporate any desired scattering probability function 
in the model. The functions included in Code 1 of this appendix are, 
the upper and the lower bounding, measured scattering probability 
functions shown in Figure 2-8. These functions are represented in 
the code by S CATER 1, and SCATER 2 respectively. Code 2 of this 
appendix is designed to handle the calculated scattering functions, 
specifically, the code includes the upper bound of the calculated 
functions shown in Figure 3-15. SCATER 3 represents this function. 
The out-puts of both codes are (1) the probability weights of each 
emerging photon, and (2) the angles of emergence. The sum of the 
probability weights for each angular range, normalized to the number 
of incident photons represents the upwelling radiances shown in 
Figures 4-1 and 4-2. 
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RADIATIVE TRANSFER COMPUTER PROGRAM 
Code 1 


PBECEDING page blank not ED.MED 
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montecarlo program with COCUMENTATION. 

TANK BOUNDARIES AND TOTAL pEPLECTIDN ARE INCLUDED. 

COMHON/BLOCKl/XMAX, YMAX,ZMAX,X,Y,Z,T,GAMA,TETA,FI.PI,DTRC,S, IS.Z-‘> 

READ DATA FROM THE T.A'^A FILF, 

PEA0( 5,25)MAXNPH,NMAXtIS 

25 FORMAT! 3 (2X. 112 ) ) 

READ! 5,3G)TETA1 ,FI1 

30 FORMAT!2!5X,F8.3)> 

READ!5,35)XMAXtYMAX,ZMAX 

35 F0RMATI3! 5X, F8.3) ) 

READ! 5,37)S 

37 F0RMAT!F8.3) 

READ!5,24)A400, A500,A600 

24 FORMAT! 3! 5X,F8.3)) 

WRITE!6,26) MAXNPH 

26 FORMAT! *0* ,• MAXI MUM NO. CF PHOTONS TO BE TPACED= *,I9) 
WRITE(6,27)NMAX 

27 FORMAT! 'O' MAXIMUM NO. CF EVENTS FOR EACH PHOTON= ’,112) 
WRITE!6,29>!S 

29 FORMAT! *C ,» INITIAL SEED FOR RANDOM NO. GENERATOR= *,112) 
WRITE!6,31)TETAI,FII 

31 “FORMAT! 'O' INITIAL TETA IN CEGREES= ',F8.3,‘ INITIAL FI INDEGR 
1EES= »,F8.3) 

WRITE 16,36) XMAX,YMAX, ZMAX 

36 =ORMAT( 'O', 'TANK DIMENSIONS IN METtRSi',' XMAX=',F8.3, 

1* YMAX=* ,F8.3, • ZMAX=',F8.3) 

WRITP(6,38) S 

38 F0RMAT!'0*, 'SCATTERING COEFFICIENT IN INVEPSED METERS= *,F8.3) 
WRITE ! 6,23) A4C0, A 50C,A600 

23 FORMAT! 'O', • absorption COEFFICIENTS AT 40C, 500, 600 NM IN INVERSE 
IMETERS:',' A400=' ,F8.3, ' A500= ' , F8 .3 , ' A600= ', F8 .3 ,////// ) 

RNW = 1 .334 


RNW is THE REFRACTION INDEX OF WATER. 

PI = 3. 141592654 
DTRC=PI/18C. 

XMAX=XMAX-*S 

YMAX=YMAX>-S 

ZMAX=ZMAX*S 

TETAI=T ETAr<D’''RC 

FII=FII-iDTRC 

NPH=1 

10 IF!NPH .GT. MAXNPH) GO TO 20CC 
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NPH IS THE NO. OF PHOTONS AT A GIVFN TIME. 

RECORD NO. OF PHOTONS TRACED AND TEST FOR END OF COMPUTATIONS. 
INITIALIZE THE COORDINATES OF THE PHOTON ENTERING THE MEDIUM. 

TETA=TETAI 

'"I='=II 

X=0. 

Y=0. 

Z=0. 000001 

DFCIDF HOW FAR PHCTCN TRAVELS BEPORE AM EVENT OCCURS. 

CALL RANDNOnS, RHCD) 

T=-ALOG(RHOD> 

GAMA=T 

T IS THE DISTANCE IN SCATTERING LENGTH UNITS PHOTON TRAVELS TO thE 
EVENT PHOTON IS AT. 

X=X+T*SIN(TETA)*CCS(FI) 

Y=Y+T*SIMTPtA)«SIN{FI) 

Z=Z+T*CQS(TETA) 

GO TO 150 
100 NPH=NPH+1 

EITHER ABSORPTION HAS OCCUR EDj OR PHOTON HAS COME OUT OF WATER, the 
PORE, START A NEW PHOTCN. 


GO TO 10 
150 CONTINUE 
KMIN=2 

IF (Z) 400,500,500 
400 XINT=X-Z*TAN(TETA)»COS(FI) 
YINT=Y-Z*TAN( TETA)*S! N( FI» 

DINT = SQRT(XINT*5f2+YlNT»«2) 
DOINT=DINT/S 

IF(D0INT .GT. 0.20)GG TO 100 
IF (RNW‘SIN( TETA ) .6T. l.C) GO TO 604 
TETAAR=ARSIN( RNW-«SIM (TEJA) ) 

1F(TETAAR .GT. 1.0)GC 100 

XINT=XINT/S 

YINT=Y!NT/S 

DINT=DD1NT 

ACT=ABS(COS(TETA) ) 

TCUT=(ABS(ZR)-ABS(Z) ) /AC"^ 

GAMA=GAMA+TCUT 

GAMA=GAMA/S 
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410 

420 

109 

838 

101 

C 

c 

c 

5999 

604 


500 

2000 

5000 


WRITE (6,410) OINT,TETAAR 

FORMAT( ///,2X,« DISTANCE FROK AXIS= * , F8 ,5 , 5X , 'POLAR ANGL E= 
WRITE (6,420) FI,XINT,YINT 

PQPMAT (' ','AZIM ANGLE= • , F8.5 , 5X , • X INT= ' , F8.5 , 5X, • YINT= 

WRITE (6, 109) GAMA 

<=0RMAT( 'O' , «GAMA = SE12.3) 

WRITE(6,888) J 

FORMATC ','N0 OF EVFNTS= • , 18 ) 

WRITE(6,101)NPH 

=ORMAT( 'O' , ' NO. OF PHOTONS TRACED = *,T8) 

CALCULAT': PHO'^ON PR38A3ILI‘^Y WEIGH*. 

CALL PHOW(PI,GAMA,DINT,A400,A500,P600) 

WRITE(6, 5999)1$ 

FOrmAT( » RANDOM NUMBER USED ’,:i2) 

GO TO 100 
KMIN=J+1 

ACT = ABS(CQS(TF-^i) ) 

TCUT=(A6S(ZR)-ABS(Z) )/ACT 

GAMA=GAMA+TCUT 

T5'^A=P1 -TETA 

FI=FI+°I 

IF(ft ,gE. 2.*PI jEZaFI-Z.^PI 

X=XINT 

Y=YINT 

Z=0. 000001 

CALL RANDNGC IS, RHOD) 

T=-AL0G(PH0D ) 

X=X+T*SIN(TETA) !<COS (FI) 

Y=Y+T=f=SIN(T?TA) -'=SIN( fI) 

Z=Z+T*CGS(TftA) 

CALL PSIVMKMIN, NMAX , J , IRTCOC ) 

IF(IRTCOD .EQ. DGO '^0 100 
!=(IRTCOD .EG. 2)G0 TO 400 
GO TO 100 
write (6,5000) IS 

format (« S'LAS"^ RAMDU'-) NUMBER U£ED = ',T12) 

‘=:T0P 

END 
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SU5R0UTINE ^SI W ( K^I N, NMAX , J ,!RT COD) 


THIS SUBROUTINE URL BE CALLED ONLY WHEN PHOTON IS STILL IN WATER 
(WHEN Z>0). 

IT DETERMINES THE COORDINATES OF THE END POINT IN THE NON-ROTATED 
SYSTEM BY =IRST ROTATING THF SYSTEM USING ANGLES TETA AND FI. 

IT GENETATES THE ROTATION MATRIX, WITH THE CONSTRAINT THAT YSTAR- 
AXIS LIES IN A PLANE PARALLEL TO THF YZ-PLANE. 

THE TCTATION MATRIX IS DESIGNATED AS A I J ( 1=1 , 3, J = 1, 3 ) . 

COMMON/ EL0CK1/XMAX,YM AX, ZMAX ,X , Y , Z , T , GAMA ,TETA , F I , PI , DTRC , S , I S , ZR 
IRTC00=0 

DO 1290 J=KMIN,NMAX 

CT=CGS{T=ta) 

C==COS(FI ) 

CT2=CT*CT 

CF2=CF*CF 

ST=SIN(teTA) 

SF=SIN( FI ) 

ST2=ST*ST 

SF2=SF*SF 

SS1=CT2+SF2»ST2 

SS=SQRT(SS1) 

SSD=1 ./SS 

A11=SQRT(1.-CF2*ST2) 
a 12=-S F»r F-^ ST2*SSD 
A13=-CT’*ST=f=CF*SSD 
A22=CT=»:SST 
A23=-SF^ST«SSD 

A31=CF*ST 

A33=CT 

A32=SF»ST 

ROTATION MATRIX HAS BEEN GENERATED. 

SCATTERING HAS OCCURFC. 

CALL ANGELS FIP,TETAP TO DISTINGUISH FROM FI, TETA 

FI®,TETAP ARE DETERMINED IN SYSTEM WITH Z-AXIS PARALLEL TO THE 

INCIDENT DIRECTTCN. 

CALL RANDNO(IS,RHCFJ 

FIP=2.*P1=*'RHQF 

CALL RANDNO{ IS,RHCT) 

CALL SCATP1(RH0T,TETA) 

T6TA=TEta*DTRC 

TETAP=TETA 

DETERMINE HOW FAR BEFORE AN EVENT OCCURS, IN THE ROTATED SYSTEM. 


77 


ORIGINAL PAGE 3S 
OR POOR QUAIZIY 



CALL RANDNO( IS, RhOD} 

T*-ALOG(RHOD) 

C 

C CALCULA'-E COORDINATES OF END POINT IN THE ROTATED SYSTEM. 

C 

XSTAR=T*SIM(TETAF)«CCS( FIP) 

YSTAR=T'*SIN{ TET4p)!rSIN( FJP) 

ZSTAR=T^CQS(TETAP) 

C 

C APPLY ROTATION MATRIX TO DETERMINE THE COORDINATES OF THE END 

C POINT IN A SYSTEM PARALLEL TO th£ ORIGINAL ONE BUT DISPLACED. 

C 

XR=A11*XSTAR+A31*ZSTAR 

YR=A12*XSTAR+A22*YSTAR+A32~ZSTAP 

ZR=A13#XSTAR+A23*YSTAR+A33*ZSTAR 

r 

C CALCULATE TETA,ANC F! IN THE PRESENT SYSTPM, WHICH IS PARALLEL 

C THE ORIGINAL ONE. 

C 

F!=ATAN( ABS( YR) /ABS( XR) ) 

IP (XR .LT* D.O)GO TO 133 
TP (YR) 333»333,633 
333 FI=2,«PI-FI 
GO TO 533 
o33 FI=FI 

GO tq 533 

133 IF (YR) 233,253t^33 

233 FI=FI+PI 
GO TO 533 
^33 FT = pi-J=i 

533 CONTINUE 
XR2=XR^XR 
YR2=YR'<=YR 
ZR2=ZR'««ZR 
Qt=XR2+YR2+ZR2 
SQOT=SQRT(DT ) 

TETA=ARCCS( ZR/SQCT) 

C 

C CALCULATE X,Y,Z OF THE END PCINT OF THF PHOTON WITH RESPECT TO 

C THE ORIGINAL AXIS* 

C 

X=X+XR 

Y=Y+YR 

X2=X*X 

Y2=Y*Y 

DIS2=X2+Y2 

XMAX2=XMAX*XMAX 

YMAX2=YMAX=4^ YMAX 
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700 

702 


704 


701 

1290 

100 

400 

500 


DIMAX2=XMAX2+YMAX2 

IFCDIS2 -GE- 0IMAX2)GQ TO 100 

Z=Z+ZR 

IF (Z) 400,400,700 
IF(ZMAX-Z) 702,702,701 
X=X-(Z-ZMAX)^TAN(TETA)*COS( FI) 

Y=Y-( 2-ZMAX)^TAN(TETA)^STN( 

ACT=AB5(CQS(teTA) ) 

TT=T-(Z-ZMAX)/ACT 

Z=ZMAX 

CALL RANCM0(IS,RH0B) 

IF{RH0B-0. 03 >704,704,100 

CHECK THREE PERCENT REFLECTION WITH UIMIFCRM ANGULAR PROBABILITY^ 

CALL RANDNQ(IS,PHCBT) 

TETA=0.5*PI^RHOBT+0.5^PI 
CALL RANDNOdS, RH08F) 

PI=^2.*PI=^RH0BF 
CALL RANDNOC IStRHGDJ 
T=-ALQG(RHOD ) 

X=X+T^SIN{TETA)^COS ('=!) 

Y=Y+T’9'SIN{TETA)=^S1N(FI) 

Z==ZMAX+T*COS (TETA) 
t=t+tt 

GAMA= GAHA+T 

CONTINUE 

IRTCOD^I 

GO TO 500 

IRTC0D=2 

RETURN 

END 
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SUBROUTINE PHPW ( PI , GAMA, DINT ,A400 ,A500, A600 ) 

THIS SUBROUTINE CALCULATES TEE PHOTON PROBABILITY WEIGHT FOR GIVEN 
WAV ELEN GTHS. 

TIR=0.0254 
TIR2=TIR*TIR 
CK=PI*TIR2 
R=0.15 

WRITE! 6, R60) R 

860 FORMAT! '0», 'SEAM RADIUS IN iMETERS= ',F8.31 
R2 = R^ P 

0INT2=DINT#DINT 

XINT=!R2-TIR2+0INT2)/!2.’!‘DI\T> 

XINT2=XINT*XINT 
YINT=SQRT! ABSIR2-XINT2) ) 

GCl=ATAN! YINT/XINT) 

GC2=ATAN!YINT/!OINT-XINT) ) 

GC3=D!-ATAN! Yl NT/ ! A 3S ! XI NT-CINT) ) ) 

AAA = GC1*R2+GC2 i=TIR2 -OINT’SYINT 
BBB=GC1*P2+GC3!?‘T1R2-0INT”'YI NT 
BIR=R+TIR 
CIR=R-T!R 

IF!DINT .GE. 0.0 .AND. DINT .LT. CIR)AREa=CK 
IF!DINT ,GE. CIR .and. DINT .LT. R)AREA=BBB 
IF! DINT .GE. R .AND- DINT .LT. BIR)AREA=AAA 
!F!DINT .GE. BIR)AREA=0. 

E500=EX°!-GAMA*A500) 

PHOTON PROBABILITY WEIGHTS ^OR 500 NM. 

PPW500= APFA'^ESOO 
WRITC!6,861 ) RPW500 

861 RQRMAT! 'O’, 'PHOTON PROB. WT. FOP 500 NM= *,F10.6) 

RETURN 

END 
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SUBROUTINE RANDNC (IX, RNUM ) 


THIS SUBROUTINE GENERATES UNIFORM PANOQW NUMBERS BETWEEN 0 AMP 1. 


IY=IX“65539 
IF(!Y) 5,e,6 

5 IY=IY+2147483647+1 

6 RNUM=IY 

RNUM=RNUM*.4656613E-9 

IX=IY 

RETURN 

END 
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subroutine SCATRKRHOTtTETA) 

THIS SUBROUTINE DETERMINES ANGLE 'THETA* FROM A GIVEN SCAT~ERING 
FUNCTION CUPPER BCUNC). 


IF(RH0T .le. 

.150)60 

TO 

1 

IFCRHO"^ .LE. 

.200)60 

TO 

2 

IF(,RH0T .le. 

.225)60 

TO 

3 

IFIRHCT .LE. 

.250,)G0 

TO 

4 

ifcrhot ,LF. 

.275)60 

TO 

5. 

IFCRHOT .LE. 

.300160 

TO 

6 

IF(R.H0T .LE. 

.320)60 

TO 

7 

IFCRHOT ..LE. 

.345)60 

TO 

8 

IFCRHCT .le. 

.360)GC 

TG 

9 

IFCRHOT .le. 

.385)60 

TO 

10^ 

IFCRHOT .LE. 

.480160 

TO 

11- 

IFCRHOT .LE. 

.550)60 

TO 

12 

IFCRHOT .LE. 

.600)60 

TO 

13 

IFCRHOT .LE. 

.655)60 

TO 

14 

ifcrhot .le., 

.6 85)60 

TO 

15 

IFCRHCT .LE. 

.715)60 

TO 

1'6 

IFCRHOT .LE. 

.730)60 

TO 

17 

IFCRHCT .,LE. 

.755)60 

■*■0 

18 

IFCRHOT .LE. 

.800)60 

TO 

1^ 

IFCRHCT .le. 

.830)60 

TO 

20 

IFCRHOT .le. 

.890)60 

TO 

21 

IFCRHCT ,LE. 

.918)GC 

TO 

22 

ifcrhqt ,le. 

.935)60 

TG 

23 

I'=C’RHOT .le. 

.945)60 

TO 

24 

IFC-RHQT .LE. 

.960)60 

TO 

25. 

IFCRHOT .le. 

.9'67')GG 

TO 

26 

IFCRHOT .LE. 

.974)60 

TO 

27 

ifcrhot ,lf. 

.981)60 

TO 

28 

IFCRHOT .LE. 

.988)60 

TO 

29 

IF-COHOT .LE. 

.994)60 

TO 

30 

ifcrhot .lb. 

1.00)60 

TO 

31 


1 TETA = C.l- 
GO TO 50 

2 TETA = aO+ (RH0T-.15>#{ .20-. 10)/ (0 .20-..l;5), 

GO TO 50 

3 TETA=.20+(RH0T-.20)vC.30-.20 )/C .225-.20) 

GO TO 50 

4 TETA=.30+(RHOT-.225)'«{.40-.,30)/( .250-. 225) 
60 ■*'0 50 

5 tpta=.40 + (RHCT- .250)»( -50-.40) /< ( . 275- . 2 50 ) 
GO' TO 50 
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TETA=.50+{RH0T-.275)*(.6C-.50) /( .300-.275) 
GO TO 50 

7 TETA = .60+( RHQT-.300)!?(.7C-.60) /(.320-.300) 
GO TO 50 

8 TFTA=.7C+( PHOT-. 320)^ (.80-. 70) / (.345-. 320) 
GO TO 50 

9 TETA=.80+{ RH0T-.345)*(.90-.80)/ ( .360-. 345) 
GO TO ?0 

10 TETA=.90+(RH0T-.360)=p(l .0-.90)/ (.385-. 360) 
GO TO 50 

11 TETA=1.0+(RH0T-.385)=i'(2.-l. )/( .480-. 385) 

GO TO 50 

12 TETA=2.0+(RH0T-.480)»(3.-2. )/( .550-. 480) 

GO TO 50 

13 TETA=3.0 + (RH0T-.550)=e<(4.-3.}/( .600-. 550) 

GO TO 50 

14 TETA=4,0+(RH0T-.600)*(5.-4. )/( .655-. 600) 

GO TO 50 

15 TEta=5.o+(RHOT-.o55)*(6.-5. )/( .685-. 655) 

GO TO 50 

16 TETA=6.0-r(RHOT-.635)r(7.-6. )/( .715-. 685) 

GO TO 50 

17 TETA=7.0+(RH0T-.715)«(8.-7. )/( .730-. 715) 

GO TO 50 

18 TETA=8.0+(RH0T-.730)*(9.-8.)/( .755-. 730) 

GO TO 50 

19 TETA=9.0 + (RH0T-.755)»-{10.-9.)/( .800-. 755) 

GO TO 50 

20 TETA=10.0+(RHOT-.800)*( 15.-10. )/(. 830-. 800) 
GO TO 50 

21 TETA=15.0+(RHDT-.83 0)>*'( 20 .- 15 . ) / ( . 890-. 830) 
GO TO 50 

22 TETA=20.0+{RHOT-.890)*(25.-2C. ) / ( .918- .890 ) 
GO TO 50 

23 TETA=25. 0+(RH0T-.91 8) *(30. -25. )/( .935-. 918) 
60 TO 50 

24 TETA=30.0+( RHOT-.935)*(35.-30. )/ ( .945-. 935) 
GO TO 50 

25 TETA=3 5.0+( RH0T-.945)=t=(40.-35.)/ (.96C-.945) 
GO TO 50 

26 TETA=40.0+( RHOT-. 9t>0 ) *( 45 .-40. ) / ( .967- .960 ) 
GO TO 50 

27 TETA=45 .0+ ( RHOT- .967 ) *( 50 .-45 . )/ ( .974-. 967) 
GO TO 50 

28 TETA=50.0+(RHQT-.974J*(60.-50.)/( .931-.974) 
GO TO 50 

29 TEtA=60.0+(RH0T-. 981)^(70.-60. )/( .988-. 981) 
GO TO 50 


OF POOB QUALEEw 
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TETA=70.0+( RH0T-,988)*( 8 0.-70. )/( .994“. 988) 
GO TO 50 

31 TeTA=80.0+(RH0T-.994)*( 180. -80 . ) /{ 1 . 00-. 994 ) 

50 RETURN 
END 
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SUBROUTINE SCATR2(RHOT,TETA ) 


THIS SUBROUTINE DETERMINES ANGLE ’THETA' FROM A GIVEN SCATTERING 
function (LOWER BOUND), 


1P(RH0'^ .Lb. 

,000)GO 

TO 

1 

IF(RHOT .LE. 

-010)60 

TO 

2 

1F(RH0T .LE. 

.01i-)G0 

TO 

3 

TF(RHO'’' .LE. 

.018)GC 

TO 

4 

IF(RHOT ,LE. 

,022)60 

TO 

5 

IFCRHOT ,LE. 

,026)GC 

TO 

6 

IF(RHOT .LE. 

.03i)GQ 

TO 

7 

I=(RHCT .LE. 

.034)GC 

TO 

8 

IF(RHOT .LE. 

.040)60 

TO 

9 

IFCRHOT .LE. 

.060)GC 

TO 

iO 

IFCRHOT .LE. 

•OSO)GO 

TO 

11 

IFCRHOT ,LE. 

.120)60 

TO 

12 

IFCRHOT .LE. 

.15060 

TO 

13 

IFCRHOT .LE. 

.175)60 

TO 

14 

IFCRHOT .LE. 

.200)60 

TO 

15 

IFCRHQT .LE. 

.220)60 

TO 

16 

IFCRHOT .LE. 

.250)GC 

TO 

17 

IFCRHOT .LE. 

.280)60 

TO 

18 

IFCRHOT .LF. 

.380)60 

TO 

19 

IFCRHOT .LE, 

.530)60 

TO 

20 

IFCRHOT .LF. 

.580)60 

TO 

21 

IFCRHOT .LE, 

.635)60 

TO 

22 

IFCRHOT .LE. 

.665)60 

TO 

23 

IFCRHOT .LE, 

.700)60 

TQ 

24 

IFCRHCT .LE. 

.740)60 

TQ 

25 

IFCRHOT .LE. 

.750)60 

TO 

26 

IFCRHOT .LE. 

.770)60 

TO 

27 

IFCRHOT .LE. 

.780)60 

TO 

28 

IFCRHOT .LE. 

.800)60 

TO 

29 

IFCRHOT ,lE. 

.833)60 

TO 

30 

IFCRHOT .LE. 

.860)60 

TO 

31 

IFCRHOT .LE. 

.885)60 

TO 

32 

IFCRHOT .LE. 

.950)60 

TO 

33 

IFCRHOT .LE. 

.970)60 

TO 

34 

IFCRHOT .LE, 

.980)60 

TO 

35 

IFCRHOT .LE. 

.990)60 

TO 

36 

IFCRHOT .LE. 

.995)60 

TO 

37 

IFCRHOT .LE. 

1.000)60 tq' 38 


1 TETA=0,2 
60 TO 50 

2 TETA=.20+(RHrT-.000)*(.30-, 20) /( .010-. 000) 
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GO TO 50 

3 TETA=.30+{RHOT-.Q10)=4=(.40-.30) / ( -014- .010) 
GO TO 50 

4 T£TA=.40+( RHQT-.014)*(.50-.40)/( -018-.014) 
60 TO 50 

5 TETA=.50+(RHOT-.018)*(.60-.£0)/ ( .022-. 018) 
GO TO 50 

6 teta=.60+(RH0T-.022)*(.70-,60)/( .026-. 022) 
GO TO 50 

7 TETA=.70+<RHDT-.026)*{.80-.70)/{ .031-. 026) 
GO TO 50 

8 TETA=.80+(RH0T-.031)^(.90-.80) /{ .034-. 031) 
GO TO 50 

9 TETA=*90 + ( OH0T-.034)*(1.0-.90)/( .040-.034) 
GO TO 50 

10 TETA=1.0+(RHOT- .040)/( . 060- .040) 

GO TO 50 

11 TETA=2.0+(RhOT-.C60)/(.090-.C60) 

GO TO 50 

12 TETA-3.0+(RHCT-.C90) / (. 120-.090) 

GO TO 50 

13 TETA=4.0+(RHCT-.120) /(. 150-.120) 

GO TO 50 

14 TETA=5.0+(RriGT-.150)/(. 175-.150) 

GO TO 50 

15 TETA=6.0+(RHCT-.175)/(.200-.175) 

GO TO 50 

16 T=TA=7.0+(RHOT-.200)/{.220-.200) 

GO TO 50 

17 teTA=8.0+(RH0T-.220)/ (.250-.220) 

GO TO 50 

18 TETA=9.0 + ( RHOT- .250 )/ ( . 280- .250 ) 

GO TO 50 

19 TETA=10,0+{RHQT-.2S0)*5./(.380-.230) 

GO TO 50 

20 TETA=15.0+(RHOT-.380)*5./( .530- .380) 

GO TO 50 

21 TETA=20 .0 + {RH0T-.530)i5./( .580-. 530) 

GO TQ 50 

22 TETA=25.0+( RHOT- . 58 0) *5 - /( . 635- . 580 ) 

GO TO 50 

23 TETA= jO. 0+( RH0T-.e3 5) «5 . / ( . 665-. 635) 

GO TO 50 

24 tetA= 33.0+( RriOT-,665)A5. /( . 700-.665) 

GO TO 50 

25 TETA=40.0+( RHOT-.700)*5./ (. 740-.700) 

GO TO 50 

26 T5TA=45.C+(RHnT-,740)*5./(.750-.740) 
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GO TO 50 

TETA=50 .0+(RHCT-.750)45./(.770-.750) 
GO TO 50 

28 TETA=55 .0+( RHOT -.770 ) *5 . / ( .780-.770 ) 
GO TO 50 

29 TETA=60.0+(RHOT-.780)*10./( .800-.780) 
GO TO 50 

30 TETA=70 .0+( RHOT-.800) *10./( -SSa-.SOO) 
GO TO 50 

31 TETA=80 .0+(RH0T-.833)*10./( .860-. 833) 
GO TO 50 

32 TETA=90.0+(RHOT-.860)*10./( .885-. 860) 
GO TO 50 

33 TETA=100,+( RH0T-.885)*10./( .950-. 885) 
GO TO 50 

34 TETA=110-+{RHOT-.950)*10./( .970-.950) 
GO TO 50 

35 TETA=120.+(RHOT-.970)*10./{ .980-.970 ) 
'go to 50 

36 TETA=130.+( RHOT-.980) *20 ./ ( .990- -980 ) 
60 TO 50 

37 TETA=150.+(RH0T-.990)*15./( .995-.990) 
GO TO 50 

38 TETA=165.+(RH0T-.995)*15./ ( 1.00-.995 ) 
50 RETURN 

END 
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RADIATIVE TRANSFER COMPUTER PROGRAM 
Code 2 


PRWEDING not pn,MED 



















A 
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o o 


C 

C 


1 '«TrCA»^Lw rtlTH JGCU'U-i\T*XT I 

'"t.'K PS AfVj TUTAL f<r,FLcCTrvi AhT IinCLUD^'. 

‘#0i ''ILis/ j:>I. OCK 1 / X MAX j Y 'lAX t Z lAX jX|Yr7^*^t»>Ar 

.CAL- a VALU 

^ AOAT/AN‘;l( 50) ,VALU( 50 ) 

'cAC )'ta r li^ SLATTr:.^ IfJO FU. 1C T i -J*'] • 


r 

c 


M ') U1,3A 

• - 1 )Af.';Ln ) ^ v^alu( i > 

il r " Hi >mT( P 10.^, J 

XFAJ JaT>. TriE MFOTUw .AT A FILL. 


i^rAol 'J, 2 MAXjJPH t >JHAX,I5 
F JP UT( J{ 2X, 112 I } 

F^AO(:>, 30) TrTAI ,ril 
33 POP 'iAl(2(pXf F8*3) ) 

-’EAO( b, 35) XMAX» Y;mX, ZMAX 
35 )R.'^AT ( 3 (5X,F8.3 ) ) 

' “ A , I ( -3 , J 7 ) s 
r J5^UT(Fd.3) 

vrVo( 5,2^ )4-'t00r A3 )3i ^\60C 
i2-+ ? (^ { 5 Xt rci. 3 ) > 

-R I T: (6t?o) ^AXN'^n 

£6 F . UT jN ’ -''iAX [MU^! jO. ,F PH'iTO^lS Tvj r^y T-^'rO^ *flj) 
w^ ITP ( 6, 2 7) 'MAX 

? 0-. '^T( » 0* , MmAMU.I iO. ,F -VFNT5 F;,- * , I li ) 

- I T t { ^ , 2 } ) I S 

>:> - AT( * OS ' IM TIAL S-E ) F )'s -.ANOC-^ nJ. GEJErv^Tv^ S I i/ ) 

^KITi= {ot -DTETA S FI I 

31 FOH'iaT{^jSM^TTI/iL TciT\ iJcGRFES- S^^*3i* ilfTI^L FI I i J- 
IEFS= SF«.3J 

-^RITF{ :»,3o)X.'IAX, Y 1AX,ZMAX 

3o f ^[</AT { »0S » T AivX .I-1*=h3iJIS I J Mr S ' aM'.a=S» 3. t 

1» Y’AX=Sf5.it» Z'UA^SPa*3) 

<KI T: (0, r d)3 

J8 FifR T( *OS ’ SLmTT:kINO Cs)rFF IC ZFjMT IN SFi.r,) 


iK I Tf (-? f 22 ) A^ )0 f AjtJJ f A600 

r''sf,,\T>( ' ,j*r'‘;\L3 j^rT J-)7'i(?jQFF-^(G‘I,eNT,Sf AT 400, 1.0J, .OJ i.” I'.V. '. 

^\L‘^7i-'S:* , • AM-).}-:* , f-S.3,» A50d=‘t ,F8.3, ’ A600= ', l-t,. :.,////// ) 

' -i'iiN= 1.3.54 


p\.-.' 13 T-i: ''cFkACTION index OF WATFP. 



90 



r-> o o o o o o \ ^ o 


Pl=3. 141592654 
^'TRC^^PI/IHO. 
a.‘!AX = XMAX S 
YMAX=YMAX' S 
/■1AX = ZMX'S 
T£T4I=T=TA1 OTrC 
FII=FII ' JTPC 
IPH=L 

1 ) 1F(<'jPH .OT. MAXMPHIGl) to 2000 

c 

iPH IS thI :-n. DF OHUIGNS AT A GIVEN TIME, 

t'EC.JrcD .u>. ■'■F PHOTJNS TRACED AND TFST FOP END Ch CmPjTaTIiu. 
I'JITIALIZF Tur Ctn^OliiATES OF THE PHOTON E-.-TFr I . , Th-- /iZOiU*. 

TA^TIITAI 
f I=FII 
X=0. 

Y=3, 

7=0. 0000 ')1 

i1 3R F'f' PHOTuN TRAVELS OPpOKt ii! EVLuT .^L'Orc>. 

:ALl kANONOU S,i?HUD) 

^=-ALjO( .\HJU) 

OAMA =T 

T IS THE 01 STANCE 1 .0 SCATTt<'I.:G LFf.GTH jMTS Pll.T./J TrsAVELS T' 
i=VLNT PmJTOv IS at. 

X=XfT'SIi'iITFTA) CJS(Fl) 

Y=Y^-^ <SI\{TET A) ‘ SIN( FI) 

2=Z+T’CuS (T5T A) 

GJ TO 150 
100 N|Ph = OPH + l 


r-ITH = P AGSoRPTIi'N mas uCCU-;tD,C’R PumTO i HA-. Cu*-*^ u JT ff 
FOFE.STAPT A i^EW PmOTON. 


SO T'J 10 
150 CONTI N0= 

KM I N= 2 

IF (/) 'tUO,500,500 
40C Xi;jT=X-Z TA’v(TFTa) cost FI) 
YINT=Y-Z- TANITETA) -SIN(FI ) 
rJINT=SO^T (Xli'IT ' 2+yiNr- ■ 2 ) 

DO I NT =D I n/S 

IF (031 IT .GT. J.2J)GJ TO 1)0 
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*+lJ 
+? J 
lo-i 

lOL 

f 

C 


500 

^JOJ 
5 000 


iF (K'W-ilfHTFTA) .GT. 1.0> GO TO G04 
TiiT AAR=AHS IN n^NU' 5 IN {T£TA ) ) 

IFtT^TAAi^ ,GT. 1.0)GO TfJ 100 

XI.NT=XINT/S 

YI.'1T = YINT/S 

::INT-0DI rjT 

/CT=A3S(Cu5(T‘=TA) ) 

TCUT = (ArvS Uh »-AOS( Z> )/ACT 
'jA'1m-=GA^ h+TCUT 
GA-^A = 0 A'iA/j) 

(>«Alo) 01 f'J T I TE T AA<' 

roU'i«Tt///,?X,'DISTAuC!i FRJ ' • , FO . 5 , -:x , 

^KITil n^O) FI,XINT,YINT 

"tV>'AT (♦ • , • ,21 1 AWGLF= ' , ?8 .5 , 5X, • a I >.T= ',Fo.-5,5a, 

iIaITC (6t 10'!?} oAMA 

FiJP-'UK *.)• ,'GA'IA = '.E 12 . 3 ) 

-P IT, (6, ,o3»J 

rilXMATC ','Ni JF £VojTS=',I 3) 
fPIT?:(ra, 101> 'IPH 

FLir« -P.K » 0» ,'M ). GF PtfOTHli TkACF^ = '.lO) 

CAICOLATC FH'jTO.J ^kGoAOIlITY vjFIGHT. 

0 1 L L PH P‘rt(PI 7lA'1AtOINTt At 00 » A3 00 i A60 U ) 

Wh ITF( o, 599V) iS 

r-UK-lATC RA ) ).) 1 aU^dlP, HSEl *,112) 

GJ ro UJ 

K'IIM=J +1 

.•CT=AS3( ;C3(TFT^I 1 
TCUT=(AoS (ZR)-A 1S(2) ) /ACT 
GANA=GA lAfTCUT 
T£Ta-PI-T5TA 

IF(FI .GC. 2. PI )‘=r=FI-2.»PI 
X=XI |T 
Y=Y I JT 
2=0.00)001 

CALL Pm'4l)iNu { 1 5 f kH.jD) 

T=-ALOo( KHiJO ) 

X = X+T"Sr.(TETA> COS(FI) 

Y=YtT"SlM(TcTAI SHKFI) 

Z=Zy-T-CC15 (TfTA) 

CALL PSrW(K(!IN,NrlAX, JjIRTCOD) 

IFIlPTCm .Eg. 1 )GlJ TO 100 
IFUPTC'jn .cO. 2JGO TO tOO 
GiJ TO 1)J 
vftITF (o,5000) IS 

rOR’lAT (’ 't'LMST •^A-'iOuM iNUMStR USED=*,I12) 


.'>L' = • ,r-,.pi 

' YI 47 = ' ,r -<.jI 
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STOP 

END 
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5Uc^0UTi:-)n PS I J (K 1IN,NMAX, J » IFTCUU) 

This SUriROUTIMc WILL 3E CALLED •'i.MLY WilEN PHOTCN ib 5 T I Ll Ii., - T I ‘ 

Z>J). 

IT DeTE-?.‘lINES TH£ CUOP:dIN AT tS OF THE ci\0 POINT IM Til*: •. i,— wT-Tf^ 
SYSTEM TY FIA5T NOTATING THE SYSTEM USING ANGLES TcTA A.,0 FI. 

IT G-NITATtS THE POTATION MATRIX , WITH THE CONSTRAINT THAT YSTA-- 
AXIS LIES IN A PLANE PARALLEL TO THE YZ-PLANE. 

THE TOTATIOh MATRIX IS DESIGNATED AS AI J { I = 1 , 3 , J= I , i J . 

DEAL'S VALU 

COMMON/ SC AOAT/ A -IGL ( 5 3 ) , V Al U ( 50 I 

common/ .-.LJCK I /X iAX,YMAX,ZMAX,X tY , Z,T ,OaiMA,TETA, f I , PI , oTr C , - , IS 
IPTC0D=3 

CO 12'J0 J=K>MI i.rMMAX 
CT = i.OSr s^TAl 
CF = COS(r i I 
CT2=CT^CT 

:f2=cf cf 

ST = 5INITF.TA1 
„F=SI-n(H ) 

3T2 = 3T''3T 

>F2?SF- 3F 

S31 = uT2-PSP2- ST2 

SS = 3 jRTISSn 

SSD=1,/5S 

1=34RT( 1.-CF2 ST?) 

4l2=-SF CF^ST? SSO 
A13=-CT' ST CF SSO 
A22 = CT* S5D 
Y23=-SF ST-SS'O 
*oi=Cr sr 

A3 3=CT 
A32^5F ST 

-sOTATIiTI ;.ATkIX has seen generated, 

SCATrcklHG HAS OCCUR EO. 

CALL ANGELS r-IP,TETAP TO DISTI-IGUISh FPOI FI,TFTA 

FIP.TCTAP A" p OETEKMI.jEO in system with Z-.iXlS PA;^ALLcL T T Tl.c 

INCIDENT OIPcCTIO-j. 

CALL PAnDimOU S, kHOF) 

FIP = 2 .'’PI . RhOF 
CALL RA^JDNUIISfRH'JT) 

Call SCATP3 (RHOT,TtTA) 

TFT i = Tl'T,j »on<c 

TETaPsTFT A 


PRECEDING^PAGE BLANK NOT FtLMEK 


94 



V ^ c ► o ri c 1 o o 


OLTc'<'‘tnE rifjW FAN BFFGRF A'J EVc JT iJCCURS,IN THE ROT-T-O SY:>Tt 


CALL RAfMutHJ( IS, kHGD) 

T=-ALJblPHOO) 

CALCULATL COLROIWATFS 0^ EM!) P-JI'MT IN THE .'LT*',TED SYSF.f’. 

XSFAR-T SI.\(T-TAP)< CGS(FIP) 

YSTA<=T ilNITcTAPI-SINIFIP) 

7STArx=r CUS(T.TAP) 

APPLY 3'UATnJj .'IrtTkIX TiJ DETFkNIWc ThE CONk JINAT=S IF Tut E''0 
PGI nT in a SYSTEH parallel to the OPI<;ifJAL CUE RUT OISPL'iCi), 

aP = a1L X3TaM- A31 'ASTAk 
Y R = 4 1 2 X 3T A,"< <• A22f Yi>TAR+A3 2 ■ Z 3 T Ar‘, 

Zk. = A 1 3* X 3 TA'^'P 425-^ YS TAP+ A3 3 ZSTAk 

CALCULAT'^ TETA.AUD FI IiN THE P^ESEHT 3 YS T F. -1 , * HICU IS PmP.'Li-EL 
THE ORIGINAL 'DUE. 

n=ATAN( ABS(Y-.)/AB3 (XR ) ) 

IF (XR .LT. J.OIGO TU 133 
IF (YR! 333,333,633 
3j 3 PI=2.'PI-ri 
jO TC 53 5 
613 FI = FI 

30 TO 533 

133 IF (YRi 233,233,A33 
233 FI=FI+PI 
GO T3 5 33 
433 FI=PI-FI 
533 CONTINUE 
XR2 = XR« X < 

YR2=YR' YR 
ZP2 = Zf<" Z? 
nT=XK2+YR2+Z\2 
SOnT = S J.'sT (DT) 

TFTA=iAPC 3S( ZR/SOOT) 

C CALCULATE X,Y,Z OF THE FND POINT O': THE PHCTuN '/ilTH PESPECT Tu 

THF jRIGINaL axis. 
r 

X=X+XR 
Y=Y+YP 
X2 = X< X 
Y2=Y' Y 
0IS2=X2+Y2 


4 




95 



O V ^ 


X''1AX2=X(1AX-^X;^AX 

YMAX2=YMAX>YMAX 

01 MAX2=XMAX2+YMAX2 

IF(0IS2 .GE. DIMAX2JG0 TO 100 

Z=Z+ZR 

IF (Z) ‘^00,400,700 
700 IF(ZMAX-ZJ70? ,702,701 
702 X=X-(Z-ZMAX) :<TAN{TETA)«C0S( FI ) 

Y=Y-( Z-Z.<1AX) ‘ TANITETA ) ''S IMF I,) 

ACT=ABS(,COS(TETA> ) 

TT=T-(Z-Z*4AX)/ACT 

Z=ZSAX 

CALL RAI'JUNOI IS 
IF(RH0'5-J.03) 704, 704,100 

CHECK TH^EE 3£RCE''iT REFLECTIOi WITH UiWIFPil'-. AHHUL^vr. Pr O'' Abl L I T Y . 

7Dt call RANDN0U$,RH03T) 

TETA=0.b PI«RH0BT+0.5->PI 
CALL RA■NO^.G( IS,RH0.3F) 

FI=2.-' PI ’‘RH06F 

CALL RANUNOIIS, RH’JDI 

T=-AL‘JG(RH0D) 

X=X + T SIi-.(TETA) « COSI FI I 
Y=Y+T‘SIN«TcT4) ' SINIFII 
Z=ZMAX+T CCSC^cTAJ 
T=T+TT 


731 

ga,^a=ga.ia 

12^90 

CCJNTINU- 

133 

tPTC30=l 


GO TO 500 

4 33 

IKTC'«0=2 

300 

•?ETUkN 
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SU«RUUri'4E: SCATR3(RH0T,T!;TA) 


C THIS SUbRCUTiNE; OETFRMI NFS ANGLE 'THcTA* FRUVi A GIVEN {AL.-’C^GY 

C CALCULATED FROM HIE THEORY) SCATTERING FUNCTION. 


REAL >i /ALU 

C OM W C N/ S C AOA T / AN G L { 5 0 ) , VA L U ( G 0 ) 

Di) iJ 1 = 1,33 

IF (kHGT .gL. VALUm .AND, ~HQT .LE. V ALU ( H-1 ) ) TO ZG 
13 CO^|TINUE 

20 TETA-AUGLI I ) + (ANoL( I+l )-ANGLt I ) ) ' (RHuT-VALU ( I) )/ ( VALJ( I+l J- 
IVA LiJ ( I ) ) 

RE TURN 
Ef^'u 
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D n ^ n 


SUo-^CUTilNE PHPW(PI ,.GAMA,OI'NTi, A400, A500,i A-60 05) 

THIS SU'iRiDUTI.HlZ .CALCULATES THE 'PHOTUw .P-RGB'ASTl ITY Wl IGHT f-j^ olVi.i 
WAVELENoThS,. 


T-IP=0,.32b4 

TIR 2 =TIR"TIA 

CK=P'I>'Tt>J 2 

=?=0.15 

WRITE(6,86))R 

A'iQ P(JRMA.T( >0’ ,’8lAH RADIUS IN MET£RS= ' ,FS.3) 

R2=k-" R 

UINT2=OI'JT'-OI’'iT 

XINT=(R2-TIP2 + 0INT2J/l 2. -UI'VT) 

XINT2=XIUT?XI.'IT 
YINT=S0PT (mDS (R2-XINT2J ) 

.GCl=ATAiUYl.iT/XINTJ 
GC2=ATAN(YINT/ (OINT-XINT) ) 
GC3=PI-ATAi,‘(YlNT/(A6S(XINT-0r>IT) } ) 

AAA = GC1 i R2+0':2’TTR2-blNT>YINT 
t«PQ=GCl R2tiJC3'TIR2-DIiMTJfYINT 
HIR^k+TIk 
CIR=R-TIR 

irfClNT .GE. C.O .AND. OINT .LT . CIP)APEA=CK 
IF<piNT ,LE. CIR .AND. DINT .LT. R)AREA=Cr.B 
IF(blMT .GC-. R .AND. DINT .LT. BIR1AR=A=AAa 
IFOI'JT ,G£. eiR)AKEA = 0. 

=5 0J=EXP(-GA--1A^ A^OO) 

V- 

C PHinCN PsJG^jBILITY weights for 500 NM, 

PPW500=APcA C500 
rtRlTF(6,361 JPPW500 

361 FORMAT! «0', •PMUTOsm PRuB. WT. FOR 500 NH= »,F10.o) 
RETURN 
END 
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SUB-^OUT UJE RANQsNOtlX, RKUn I 
r 
C 

THIS SUBROUTINE GENERATES UMFOR»'1 RANDOM UU 1dER:> dETwEEU 0 hi‘>‘ 
C 
r 

lY=Ix^ 63539 
IF(IY) 5ro,6 
5 iY=I Y+^1474<i3647fI 
o RNUM=IY 

-^.46 50 6135-9 

IX = IY 
Kt TURN 
rNO 
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FILE 

04; SCATR3 DATA A 


Oi J 

0.2000 

0. 4000 
0.6000 
0.6000 
1. 003 J 
1.2000 
1 .4000 
1.6333 

1. aooo 
2.0000 

10.0000 
18 .0000 
26. 330 0 
34 .0300 
42.0000 
6 3.333 ) 

58.0000 

66.0000 

74.0000 

52.0000 
136.333 3 

114.0000 
122 .0000 
13 3. 333 3 

133.0000 

140.0000 

154.0000 

162.0000 
170.0003 
178.2000 
178 .6000 
179.2033 
180.0000 


J.OOOOOOD'^OO 
0.253470D<-00 
0.577503D+00 
0.6 73157D + 00 
0. 72H562DfOO 
3.7516a60+)3 
0. 76982 9D+00 
0 .7328220+00 
3.7928120+33 
0. 8005860+00 
0.6068250+00 
0.9301360+00 
0.9493300+00 
3.9633730+30 
0.9730730+00 
0.980206D+00 
3.9850^40+33 
0.9384060+00 
0.990774D+00 
0. 9925100+00 
0.9938340+00 
3.9966930+33 
0. 9973770+00 
0.99796013+00 
3.903464D+33 
0. 9988980+00 
0.9992650+00 
3.99955 30+00 
0. 9997300+00 
3.9999270+30 
0.999997D+00 
0.9999980+00 
3.9999990+33 
0. lOOOOOD+Ol 


100 



FILE 



05 «F or ON 

DATA 

A 


12 

100 

53479 

9.700 

0. 


1.217 

12.000 

1.217 

2 . c>0 0 

4.300 

3.200 

2,400 
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file: mfotqn exec 


A 


GL TXTLIB FORTMQOl 
FI 04 DISK SCATR3 'DATA -A1 
FI 05 DISK MFOTON DATA A1 
FI 06 PRINTER 
LOAD MFQTON 


START 



APPENDIX B 


MEASUREMENT OF SCATTERING FUNCTIONS 


B.l Scattering Functions 

Scattering is an inherent property of the water which is useful 
as an optical parameter. Detailed knowledge of the scattering 
functions, in fact, can yield information about the particle size 
distribution and the composition. 

The scattering function o(6) is defined by the relation 

o(e) = (meter ^ Str (1) 


where dJ(0) is an element of radiant intensity scattered in the 
direction 0 from the incident beam by the volume element dV* H is 
the irradiance received by the sample volume. 

B. 2 Determination of Volume Scattering Function 

Both the sample volume and the small solid angle, within which 
the radiant intensity is measured, are determined by the optical 
geometry of the instrument used. The instruments usually use a 
finite sample volume and collect the energy scattered at angle 0 
over some solid angle. The equation (1) is then written as 




J(e) 

H.A.A 


P(6) 

P(0).S^.Jt 
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where 


P(0) = the total light flux entering the savaple volume 

P(0) = the Ixght flux entering a small solid angle about 
the angle 0 at which the measurement is made 

9, = the solid angle over which the measurement of P(9) 

is made 

£ ~ length of sample volume 

A = the projected area of the sample volume V, as seen in 
the direction of P(0) 

It is necessary to know either P(0) or P(0) in absolute terms. 

The scattering instruments allow P(0)/P(O) to be computed. The 
length, and the solid angle, 9 are determined by the geometry of 
the instruments. 

l^en a scattering measurement is made using a finite volume of 
water, an unavoidable error is caused by absorption in the sample 
volume. If the instrument used had a sample path length that is 
small relative to the attenuation length of the water, this error 
is small and is less than the instrument errors. If the measurement 
is made using a path length that is not small relative to the 
attenuation length, the results have to be corrected. One such 
correction applied can be referred in Reference 10. 

B. 3 Scatterance Meters 

The scattering quantities have exact mathematical definitions 
which dictate the design of the meters to be used. In principle, 
measurements of scatterance involve Irradiation of sample volume by 
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a beam of Ixght and recording of the light scattered by the volume 
through various angles. 

Several types of scatterance meters have been developed. 

Typical types are: Fixed angle, Free angle, and Integrating meters. 

One ground of subdivision is to distinguish in-vitro and in-situ 
meters. 

It is not our intent to describe in detail various scattering 
meters used by researchers in this area, however, a brief discussion 
may be warranted regarding the differences between general type and 
small angle scattering meters. Typical scattering meters used by 
Scripps Institution of Oceanography are briefed below. 

B.3.1 General Angle Scattering Meter 

Its purpose is to determine volume scattering function between 
the limits of 10^ < 0 < 170^. It has a projector which rotates 
about the sample volume from 0=0^ through 0 = 180^. The measure- 
ment at 0 = 0^ indicates total power in the projected beam, while 
the measurement at 180*^ records the background ambient light level. 
The rest of the readings (10° < 0 ^ 170°) measure scattered light. 

The output of this instrument contains analog voltage signals 
representing depth, scan angle position and the photometer signal. 

A continuous trace of photometer signal versus depth at any fixed 
angle between 10 and 170 degrees and a continuous trace of photometer 
signal versus scattering angle at a fixed depth are the two methods 
of data collection using general angle scattering meters. 
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original 
OF POOR Q,UALXry 



B.3.2 Small Angle Scattering Meter 

Small angle scattering meter (with, which the results used in 

Reference 1 were measured) xs essentially that which was modified 

(3) 

and used by Morrison. Main problems in low angle scattering 
meters are; scattering within the instrument and limitations in 
defining the limits of solid angle of the measurements. 

In the low angle scattering meters used in Reference 10, and 
attempt is made to reduce the instrument's over internal scattering. 
It is still significant relative to the small angle forward scat- 
tering of clear waters. 

The instrument has a projector which having a small point source 
of light, produces a beam of highly collimated light. After tra- 
versing the sample path, the light enters a long focal length lens 
in the receiver and an image of the point source is formed at its 
focal length. The light which traverses the water and is neither 
absorbed nor scattered falls within this small image. Light which 
is scattered arrives at the image plane displaced from the axis at 
a distance proportional to the angle through which it has been 
scattered and is the focal length of the receiver lens. The scat- 
tered light is allowed to pass through four field stops before 
reaching the detector. The first field stop is a small hole and the 
other three field stops are annulus. The inner and outer radii of 
the annulus determine the angular interval over which the scattered 
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light is accepted. The solid angle, Q, in equation (2) is liniited 


by the angles 9^, 6^ imposed by the annulus field stops and is 

2 2 

calculated, from fi = 2ir (02 - 0^), where 0^, and are in radians. 

The value computed for the volume scattering function a (6) is 
an average value for a(0) between the angular limits, 9^^, and 02 > 
of the solid angle. 


ORIGINAL PAGE IS 
OF POOR QUAI^ 
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APPENDIX C 


LISTINGS FOR POLYMIE AND DBMIE ROUTINES USED TO 
CALCULATE THE VOLUME SCATTERING FUNCTIONS 


PREX3BDIN6 PAGE BLANK NOT FUMED 



o r- r-j o o o r> c > o r> o o o o o o o o o o o o n o o 


MAIN 


MAIN PROG POL YM I El VECTOR) 

THE FOLLOWING, DOUBLE PRECISION INPUTS ARE REQUIRED: 

rfk=rfal part of refractive index 

RF 1= IMAGINARY PART OF REFRACTIVE INDEX 
RADU=UPPER8DUND ON RADIUS (MICRONS ) 
wAV£=WAVELENGTH IN MICRONS 
A{ I)=PARAMETERS for DISTRIBUTION ONE 
S( I )=PARAM£TERS FOR DISTIBUTIDN TWO 

THFTD( J)=VECTOR OF ANGLES FROM 0-90 (COMPLIMENTS ARE ALSO CALC) 
OTHER INPUTS ARE: 

JX=NUMBER OF ANGLES FROM 0-90 
NPAD=NUMBER OF RADII BETWEEN 0-RADU 
NPARA=NUM8ER UF PARAMETERS IN DISTRIBUTION ONE 
NPARA2=NUMBER OF PARAMETERS IN DISTRIBUTION TWO 
TW0=L0GICAL VARIABLE TO ENABLE THE USE OF TWO DISTRIBUTIONS 


TwO FUNCTION SUBPROGRAMS DIST(RAD,AJ AND 0IST2(RA0,B) ARE REQUIRED 
IN ADDITION TO POBMIE SUBROUTINE 


TWO DATA SETS (6 AND 8) ARE USED FOR OUTPUT; NORMALLY 6=PRINTER 
AND 8=TAPE 


10 Fi)RMAT{3015.5) 

11 F0RMAT(2015.5) 

12 rORMAT(D15.5) 

13 F0RMAT(D15.5, 15) 

14 FORMAT (21 5) 

16 F0RMAT(L5) 

16 FORMATdS) 

17 F0RMATIO15.5) 

20 FORMAT(lHl) 

25 F0RMAT(//T10, 'ELEMENTS OF THE TRANSFORMATION MATRIX FOR A SPHERE 
IWITH SIZE PARAMETER = »,F15.5) 

30 F0RMAr(//T10t'REFRACTIVE INDEX. REAL= » ,015 .5, T60, • IMAGINARY* , 015 
1.5,//) 

35 F0RMAT{T3, 'ANGLE' ,T17,' SIGMAl* ,T3 I , • SI GMA2 ' , T46 , • S IGMA 3' ,T61, 'SIGN 
1A4* ,T76, ' INTENSITY* , T9 1, ' POLARIZATION'// ) 

40 F0RMAT(F10.4, 5E15.6,F15.4 ) 

. ,V5 -Fj3RHAT-(/7T;l'J',.'>,'EFF,rCT.ENCY -FACTOR FOR EXT INCTIQN* ,E 15.6) 

‘ 50 FORMAT! />f 10,' EFFICIENCY FACTOR FOR SCATTER ING* ,E15 .6 ) 

55 F0RMAT(//T10, * EFFICIENCY FACTOR FOR ABSORPT ION' ,E 15.6) 

60 F0RMAT(//T10, ' ASYMMETRY FACTOR' , E15 .6// I 
70 F0RMAT(//T10, • TOTAL TIME FOR THIS CASE IN SECONDS* ',F15.3//) 
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2 FORMAT{//T10, 'PROBABILITY FOR THIS SIZE PARAMETER = *,015.5,//) 

3 F0RHAT(//T1J, 'NORMALIZATION FACTOR FOR THIS SET OF SIZE PARA=', 
ini5. 5,//) 

REAL- 8 RFR,RFI,X,»JEXT,gSCAT,QA 8 S,THET 0 { 100 ) ,PQEXT, PQ SCAT ,PQABS 
dO F0RMAT(//T10,' SCATTERING CROSS SECTION' , E15. 6 > 

REAL'8 ELTRilX{4, 130, 2 ) , AL AM, CON, CTBRQS , AVCSTH , PELTMXI 4, 100 ,2 ) 
REAL ‘‘8 PAVCTH,THE( 100) ,PBSCAT 
REAL'^4 AIN(100,2 ),POLRUOO,2) 

PEAL •‘4 PAIN( 100,2) ,PPOLR(lOO,2) 

REAL 4 PAI(100,2),PPOL(100,2) 

REAL -8 PR0B2,PN0RH2 

REAL 8 PQEX,PQSCA,PQAB,PBSCA,PAVCT,PELTH<4,100,2) 
kEAL 8 RADU, ORAD, WAVE, GAMMA, A(20) ,PR0B,B(20) 

LOGICAL WRNfTWO 
WRN=, FALSE. 

C0N=3. 141 5926535 897932D+0 
INTEGER NPARA,NPARA2 
90 READ (5,10) RFR,RFI,WAVE 
READ (5,14) JX,NPARA 
READ (5,12) { THETDd ) ,1=1 ,JX) 

READ (5,13) RAOUtNRAD 
1 FORMAT (015.5) 

DO 5 I=1,NPARA 
READ (5,1) A(I) 

5 CONTINUE 


PEAD(5,15) two 
DO 95 1=1, JX 
95 THE( I )=THETD( I ) 

IF (TWO) GO TO 61 
GO TO 62 

61 REA0(5,16) NPARA2 
DO 62 I=1,NPARA2 
READ(5,17) B(I) 

62 CONTINUE 
PuEXT=O.ODO 
POEX=O.ODO 
PySCA=0.000 
PQAB=O.ODO 
P8 SCA=O.ODO 
PAVCT=O.ODO 
PQ SCAT=O.ODO 
PQABS=0,0D0 
PB SCAT=0,ODO 
PAVCTH=0.000 
ORAD=RADU/NRAO 
DO 1000 J=1,JX 
DO 1000 K=l,2 
DO 999 1=1,4 
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PELTMX(I,J,KJ=0.000 
PELTM(I,J,K)=0.0D0 
CONTINUE 
PAIN(J,K)=0.0D0 
PAH J,K)=O.ODO 
PPOL(J,K)=O.ODO 
PPOLR(JtK) =0.000 
1000 CONTINUE 
PAD=0.0 
,PNGRM2=0.000 
IF (TWO) GO TO 91 
P.N0RM2=1 .ODO 

91 CONTINUE 
Pi'iORM=0. JDO 
TI ME 1=0.0 
no 3030 L=1,NRA0 
RAO=RAO+-ORAD 
DO 100 J=1,JX 
100 ThETD(J)=THE( J) 

X=2.0D0 CON- RAO/WAVE 
PROB=OIST(RAD,A) 

IF ITWO) GO TO 63 
PR0B2=0.0D0 
GO TL) 64 

63 PkCB2=0IST2(RAD,B) 

64 CALL SETCLK 

CALL POBfllE I X jRFRtRFI ,THETD, JX,QEXT,QSCAT,CT8RQS,ELTRMX,WRNI 
CALL READCL(TIME) 

If (WRN) GO TO 1001 
PNURM=PNORM+PROB 
PNORM2=PNORM2+PR082 
TI ME1=TIME1+TIME 
.JA6S=QEXT-QSCAT 
AVCSTH=CTdRQS/OSCAT 
DO 150 K=l,2 
DO 150 J=1,JX 

AI N{ JtK)= ELTRMXI I, J,KM-ELTRMX{2,J,KI 
' > POL-RIJ,KI= (ELTRMX(2,J,K»-ELTRMX(1,J,K))/AIN{J,KJ 

Allir,rjfK}= ,5'^AIN(J,K) 

PAHJI J,K)=AIN( J,K) ipROB+PAIN (J,K) 

PA I (J,K)=AINI J,K)'-PR0B2+PAI ( J,K) 

PPOLI J,K) = POLR( J,KKPRQB2<-PP0L(J,K) 
PPOLR(J,K)=PPOLR(JfKH-POLR(JtK) •'PROB 
150 CONTINUE 

DO 2000 1=1,4 
DO 2000 J=1,JX 
DO 2 000 K=l,2 

PELTMXI I , J,K) =PELTMX(I , J,K)+ELTRMX(I,J,K)'=PR0B 
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PELTM( I, J,K) = PELTM( I , J ,K ) -i-ELTRMXI I , J , K ) » PR0B2 
2000 CONTINUE 

WRITE(6, 20) 

1«RIT£(6,25) X 
wRITE(6,30) RFR.RFI 
WRITE(6,35) 

WRITE (6, 40) { (THETD(J)»IELTRMX(I , J , 1 ) , I =1 1 4) t AI N IJ 1 1 ) » PO LR { J , 1 ) ) f 
1J=1, JX) 

WRITE (8,40) ( (TH£TO( J ) , (ELTRMX( I , J , 1 ) , 1= 1, 4) , AIN { J , 1 ) , POLR ( J , 1 ) ) , 
IJ=1, JX) 

DO 2J0 J=1,JX 

THETO(J)= 180.0D0-THETD( J) 

2JJ CONTINUE 
JMX=JX-l 
00 210 J=i,JMX 
JJ=JX-J 

WRITE (t>,40) (THETDI JJ ), ( ELTRMX ( I , J J , 2 ) , 1= 1, 4 ) , AIN ( J J , 2 ) , POLR ( J J , 2 ) > 
L WRITE{3,40)( THETD( JJ) , (ELTRMXd ,JJ,2),I=1,4),AIN(JJ,2) ,POLR( JJ,2) ) 

21 0 CONTINUE 

WRITE(o,45) QFXT 
WRITE(6,50) QSCAT 
WRITE(6,55) QA8S 
WRITE(o,60) AVCSTH 
WRITE (6, 2) PROB 
wRITE(6,2) PR0B2 
WRIT£(6, 20) 

«RIT£(6,70) TIME 

PQS C AT- P QSCAT +QSCAT-' PROB 

PQSCa=PQSCA+QSCAT^PR0B2 

PQEX=PQfcX+QEXT»PR0B2 

PuAB = PQAtl+-QABS» PR0B2 

PBSCA=P6SCA+QSCAT ■ (RAD-» ‘-2 )-‘PR0B2 

PA VCT=PA VCT+AVCSTH* PR0B2 

PQEXT=PUEXT+QEXTi'PROB 

PQABS=PyA8S+aABS’'PR06 

PBSCAT=PBSCAT+QSCAT» (RAD^-«2)“'PR0B 

pavcth=pavcth+avcsth*prob 

1001 WRN= .FALSE. 

3000 CONTINUE 

DO 4000 J=1,JX 
DO 4000 K=l,2 
00 4001 1=1,4 

PELTMXn ,J,K) =PELTMX(I , J,K)/PNORM 
PELTMd, J,K)=PELTM( I , J , K ) /PN0RM2 

4001 CONTINUE 

PAIM( J,K)=PAIN(J,K)/PNORM 
PA I (J,K)=PAI(J,K)/PN0RM2 
PPOL U,K.)=PP0L{J,K)/PN0RM2 
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PPOLRCJ ,K)=PPOLR( Jf K)/PNORM 
CONTINUe 

C END FILE 8 

PUSCAT=PQSCAT/PNORM 
PyEXT=PQEXT/PNORM 
PQABS=PyABS/PNORM 
PB5CAT=P8SCAT -CON/PNORM 

pavcth=pavcth/pnorm 

PQSCA=PQSCA/PNORM2 
?QlX=PQEX/PN 0RM2 
PQAB=PyAB/PNORM2 
PP SCA = P»SCA* CON/ PN0RM2 
PAVCT=PAVCT/PNORM2 
DO 6000 J=1,JX 
6000 THE1D(J J=THE( J) 
wRITE(6,20) 

65 F0RMAT<//T10, » ELEMENTS OF TRANFORMATI ON MATRIX FOR POLYDISPERSION' 

1 ,//) 

WRIT£{6,65) 

WRITE(6,30) RFR,RFI 
WRITE(6,35) 

WRITE (6,40) ( (THETDI J) , ( PELTMX 1 1 , J » 1 ) r 1= 1 1 4 ) , PAI N ( J> 1) ,PPOLR (J, 1) 
1), J=1,JX) 

C wRIT£{8t40) UTHETDI J) ,(PELTMX(I ,J»l J = ,PAIN(J, l»,PPQLR (0,1) 

C 1),J=1,JX) 

DO 5000 J=1,JX 

THETO( J)=180. JOO-THETD(J) 

5000 CONTINUE 
JMX=JX-1 

DO 5001 J-1»JMX 
JJ=JX-J 

yR ITE(6,40» ( THETO{ JJ) ,( PELTMX (I , JJ , 2 ) ,I =1 ,4 > , PAIN(J J,2) ,PPOLR 
1 {JJ,2) } 

C ^^RITE^8,40) (TH£TD( JJ) , ( PELTMX ( 1 , J J , 2 ) , I = 1 , 4 ) , PA IN ( J J , 2 > ,PPOLR 

C 1(JJ,2)) 

5001 CONTINUE 

C END FILE 8 

WRITE(6,45) PQEXT 
WRITE(6,50) PUSCAT 
WRITE(6,55) PyABS 
WRITE(6,80) PBSCAT 
WRITE(6,60) PAVCTH 
WRITE(6,3) PNORM 
WRITE(6,70) TIMEl 
WRITE(6,20) 

DO 5010 J=lrJX 
5010 THETD( J) =THE{ J) 

IF (TWO) GOTO 5102 
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GO TO 5003 
5002 WRITE(6,20) 

hRITE(6,65) 

WRITE (6, 30 » RFR,RFI 
WRITE{6,35) 

WRITE (6, 40) ( (TrlETD(J»,<PELTMn,J,l) ♦ I =li 4) , PAI (Jf 1)| , PPQL IJ t U 
1) , J=1,JX) 

C WRITE (8, 40) { (THETOU ) , ( P ELTM (I » J 1 1 ) » 1= 1 ,4) , PAU J , 1) ,PPOL(Jf II 

C 1).J=1.JX) 

00 5004 J = 1tJX 
THETDIJ) =180. ODO-THETO( J) 

5 0 0*^ CONTINUE 
JMX=JX-1 
DO 5005 J=1,JHX 
JJ=JX-J 

WRITE (6, 40) (ThETDUJ), (P ELTM <1 1 J J f 2 ) , I= 1, 4 ) , PA I ( J J, 2)fPP0L 
1(JJ,2)) 

C WRITb{8,40) (THETD{JJ),(PELTHa,J0«2),I = l,4),PAI(JJ,2),PP0L 

C 1UJ,2)I 

5005 CONTINUE 
C END FILE 8 

WRITE(6,45) POEX 
WKITE(6,50) PQSCA 
WRIT={6,55) PQAB 
WRITE(6,80) PBSCA 
WRITE(6t60) PAVCT 
WRITE (6,3) PN0RM2 
WRITE(6,70) TIMEl 
WRITE(6,20) 

5003 STOP 
END 
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PD8HIE 


SUBRCUTINE PDBMIt ( X ,RFR, RF I , THETDt JX ,QEXT , QSCAT ,CTBRQSt ELTRMX, KRN 
1 ) 

-RADIATION SCATTERED BY A SPHERE. THIS SUBROUTINE CARRIES OUT ALL 
SUBROUTINE FOR COMPUTING THE PARAMETERS OF THE ELECTROMAGNETIC 
COMPUTATIONS IN SINGLE PRECISION ARITHMETIC. 

THIS SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE OF 
OOWNWARO RECURRENCE RELATIONSHIP. 

X 0 SIZE PARAMETER OF THE SPHERE,! 2 ^ P I ‘ RADIUS OF THE SPHERE)/ 
.WAVELENGTH OF THE INCIDENT RADIATION). 

RFD REFRACTIVE INDEX OF THE MATERIAL OF THE SPHERE. COMPLEX 
QUANTITY. .FOR MO {RFR - I * RFI ) 

THETD(J)3 ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT 
AND THE SCATTERED RADIATION. THETD(J) IS OR = 90.0. 

IF THETD(J) SHOULD HAPPEN TO BE GREATER THAN 90.0, ENTER WITH 
jUPPLEMFNTAkY VALUED SEE COMMENTS BELOW ON ELTRHX.. 

JXO TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATION AREREQUIROE. 

JX SHOULD NOT EXCEED 200 UNLESS THE DIMENSIONS STATEMENTS 
ARE APPROPRIATELY MODIFIED. 

MAIN PROGRAM SHOULD ALSO HAVE REAL THETD { 200) , ELTRMX ( 4 ,2 00 ,2 ) . 
DEFINITIONS FOR THE FOLLOWING SYMBOLS CAN BE FOUND IN • LIGHT 
SCATTERING BY SMALL PARTICLES, H. C. VAN OE HULST, JOHN WILEY + 
SONS, INC., NEW YORK, 1957 '. 

0EXT82 cFFIECIENCY FACTOR FOR EXTINCTION, VAN D£ HULST, P.14 + 127 
USCAT82 FFFIECIENCY FACTOR FOR SCATTERING, VAN DE HULST, P.14 + 127. 
CT6ROS3 AVERAGEICOS IN6 THETA) * 0SCAT,VAN DE HULST, P. 128. 
tiLTRMXd ,J,K> 0 ELEMENTS OF .THE TRANSFORMATION MATRIX F,VAN DE HUL 
iT, P.34,45 + 125. I = 10 ELEMENT M SUB 2.. I = 23ELEHENT M SUB 1.. 

I = 30 ELEMENT S SUB 21.. I = 40 ELEMENT 0 SUB 21... 

ELTRMXII , J,l) REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR 
THE ANGLE THETD(J).. ELTRMX(I,J,2) REPRESENTS THE ITH ELEMENT 
I'F THE MATRIX FOR THE ANGLE 180.0 ~ THETD(J) ,, 

FOPMATIlOX* THE VALUE OF THE SCATTERING ANGLE IS GREATER THAN 90.0 
$ DEGREES. IT IS • ,EL5,4) 

PDRMAT(//1JX' PLEASE READ THE COMMENTS*//) 

F0RMAT(//10X‘ THE VALUE OF THE ARGUMENT JX IS GREATER THE 100*) 
F0RMAT(//1JX' THE UPPER LIMIT FOR ACAP IS NOT ENOUGH. SUGGEST GET 
IDETAILED OUTPUT AND MODIFY SUBROUTINE*//) 


REAL'S X , RX , R FR , RF I , y EXT , QSCAT ,T(5),TA(4),TB(2),TC(2),TD(2),TE!2) 

2 CTBRyS ,FLTRMX( 4, 100,2 ), PI (3,100) ,TAU(3,100) , 

3 CSTHT(IOO) , SI2THT(100 ), THETD! 100 ) 


COMPLEX' 16 RF ,RRF, RRFX, WM L,FNA , FNB, TC 1 , TC2, WFN ( 2 ) , AC AP (8 000) , 
2 FNAPtFNBP 


9 


LOGICAL WRN 

FORMAT!//TiO, ’WARNING, ACCURACY NUT ACHIEVED*//) 
TA!U0.REAL PART OF WFN(l)., TA(2)0 IMAGINARY PART 
T4(3)0 REAL PART.,0F- WFN (2 > . . TA(4)0 IMAGINARY PART 


TB( 1 )0*kE‘AL 
TC(1)0 REAL 


rart 

PART 




,0'F'FNA...TB(2)0 

0FtFNB...TC(2)0 


IMAGINARY 

IMAGINARY 


PART 

PART 


OF 

OF 


OF WFN(l).. 
OF WFN!2).. 

FNA.. . 

FNB ... 


T 
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C TD(1)0 REAL PART OF FNAP.. TD(2) rMAGINARY PART OF FNAP«,. 

C TCCDO REAL PART OF FNBP... TE(2)0 IMAGINARY PART OF FNDP... 

C FNAP + FN8P ARE THE PRECEDING VALUES OF FNA + FNB RESPECTIVELY. 

EQUIVALENCE <WFNU)f TA ( 1 )) , (FNA, TBU)), (FNB, TC(U) 
equivalence (FNAP, TD(U), (FNBP, TEdH 
IF ( JX .LT, 101 1 GO TO 20 
'iRITE (6, 7) 

XRITE16, 6) 

STOP 1 

20 RF=OCMPLX(RFR,-RFI) 

RRF = l.ODJ/RF 

RX = 1.00 0/X 

RRFX = PRF - RX 

T(1)=(X' 2» (RFR->'.'2 + RFI=f'^2 ) 

r(l)=DSyKT(T(l)) 

NMXl = 1. 1000 T( 1) 

IF (NMXl .LT. 7999) GO TO 21 
v^RITE(6, 8) 

STOP 2 

21 NMX2 = T(l) 

IF (NMXl .GT. 150) GO TO 22 
NMXl - 15 0 
NMX2 = 135 

22 ACAP(NMX1 1 ) = ( 0.000, O.ODO ) 

DO 23 N = 1, NMXl 
NN = NMXl - N + 1 

ACAP(NN) = (NN-H) ' RRFX - 1 .000/ ( (NN + 1 )»RRFX + ACAP(NN+D) 

23 CONTINUE 

Dil 30 J = 1, JX 

IF { THETD(J) .LT. 3.0DJ ) THETD(J) = DABS( THETOi J ) ) 

IF ( THETO(J) .GT. O.ODO ) GO TO 24 
CSTHT(J) = 1.300 
SI2THT(J) = O.ODO 
GO TC 33 

24 IF ( THtTD(J) .GE. 90.0D0 ) GO TO 25 

T(l) = { 3.141592b535697932D<-3 THETD{ J ) )/ 180.D0 
CSTHT(J) = DC0S(T(1) ) 

SI2THT(J) = 1.3D0 - CSTHT(J)-»“2 
GO TO 30 

2p IF ( THETD(J) .GT. 90.000 ) GU TO 28 
CSTHTU ) = O.ODO 
SI2THT{ J) = 1 .000 
GO TO 30 

23 WRITE (6, 5) THETD(J) 

WRITE(6,6) 

STOP 3 

30 CONTINUE 

DO 35 J = 1, JX 
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300 

33 


6 0 


6 5 


70 


PI (If J) = 0.000 
Pl(2,J) = l.ODO 
TAU(lfJ) - 0.000 
TAU (2, J) = CSTHT(J) 
CuNTINUt 


T( 1 ) = DCOS(X) 

T(2) = f)SIN(X> 
wfll=DCMPLXm 1) i-T(2) J 
wFN( 1)=DCPPLX(T(2) iT(in 
WF.N(2) * RX WFN(l) - W'U 
TC 1 = ACAPd J •*- RRF + RX 
TC2 = AGAPd) ' RF + RX 

FNA = {TCI ■ TA(3) - TA ( I J) / (TCI * WF!M(2) - WFN{ IN . 

FNB = (TC2 » TA(3I - TA(IJ> / {TC2 WFN{2) - WFN(l)) 

?NAP = FNA 

FNOP = FNB 

T( 1} = 1. 5000 

TB ( 1 > = T d J - T B d ) 

T3(2J = T(l) - TB(2) 

TCU) = Td) ' TCdJ 

TC(2) = T(l) TC«2) 


Od 6J J =1,JX 
tLTRHXI 1, Jfl) 
FLTRMX(2, Jfl) 
FLTRMXOf Jf 1) 
ZLTRMX(4t Jfl) 
rLTRMXI If Jf 2) 
tLTkilX(2f Jf2) 
ELTRi-IXtif J,2» 
FLTRMX(4f Jf2) 


T3(l) PI(2,J) <■ TCd) TAU(2,J) 

TS{2) • PI(2fJ) +TC(2) ^ TAU(2fJ) 
TCm 'PI(2fJ) + TBdJ TAU(2fJ) 

TC(2) 5 PU2,J) + TB(2) < TAU(2,J) 

T8(U ' PI(2fJ) - TC(U * TAU(2fJ) 

TB(2) -- PI(2,J) - TC(2) * TAU(2,J) 

TC(U - PI(2,J) -TS(1) ’ TAU(2tJ) 
TC(2> - PI(2,JJ - TB(2) =• TAU(2fJ) 


CONTINUE 

'JEXT = 2.0DO ■' { TBd) + TCd)) 

OSCAT =(TBd)-''2 + T8(2)^*2 + TCt.lJ»--2 + TC ( 2 )"' 2 ) /O . 7500 
CTBRQS = 0.000 
N = 2 


T( 1) = 2 - N - I 
T(2) = N - 1 
T(3) = 2 -‘N + 1 
D070J=1 f JX 

PI(3,J) = (T{l)-PI(2fJ) -CSTHT( J)-N''PI( If J ) )/T( 2) 

TaU(3, J)=CSTHT{J)^ (PI (3,J )-PI (1, J ) )-T ( 1 I 2THT ( J ) - P I { 2f J ) +TAU( If J 


1 ) 

CONTINUE 
WMl = WFNd) 
wFN(l) = WFN(2) 

VfFN(2) = T(U » RX ' WFNd) - WMl 
TCI = ACAP(N) ^ RRF + N RX 
TC2 = ACAP(N) ‘ RF + N • RX 
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F.mA = (TCI -TAOJ - TA(D) / (TCI '' WFN(2) - WFN(l)) 

FH6 = (TC2 * TA(3) - TA(1)) / (TC2 < WFN(2) - WFN(1)> 

T(5} = N 

T(4> = Td) / (T(5> T(2)J 

T(2) = (T(2)MT(5) * 1.0D0))/T(5) 

CTORQS = CTBROS + T(2) • (TOm TB(1) + T0{2) > TB ( 2 ) + T£ ( 1 ) - 
iTC(U + TE(2) ' TC(2)J + T(4) 5 (TO(l) » TE(1) + TD(2) TE(2H 
JEXT = gEXT + T(3) - (TB(1) + TC(D) 

T{4) = TB(l)k 2 + T8(2)-*"2 <■ TC(LJ»-2 + TC(2)>'''2 
OSCAT = JSCAT + T(3) T{4) 

T(2) = N ’ (N + 1) 

T( 1) = T(3) / T(2) 

K = (N / 2) ' 2 

00 8) J = 1, JX 

FLTRMXd, J,1 »=ELTkMX( 1 , J ,1 ) +T d ) -- { TB d > •" P I ( 3 , J » +TC d d TAU ( 3 , J ) > 
£LTRMX(2, J,1 )=tLTRMX(2, J,l)+Td )--'(TB(2d PI(3,v))+TC(2) •'TAUO, J) ) 
ELTRMX(3, J,l) =£LTRMX(3, J, Udd d(TCd I ’PI( 3, J J+TBd > -TAU(3, J) ) 
LLTRRX(4, J,1)=ELTRMX(4, J, 1>+T(U--{TC(2)''PI{3,J)+T8(2) --TAU(3, J) ) 
IF(K.Eg-N) GO TO 75 

FLTRMXC 1, J ,2>=ELTRMX( 1 , J , 2) +T d ) M TB ( U • PI ( 3 , J ) -TC d ) ' TAU{3, J) ) 
-LTRMX(2,J,2dELTRMX(2, J,2J+Tdd (TB(2) > P I ( 3, J J -TC ( 2 ) - TA U { 3, J M 
FLTRMX(3, J,2)=ELTRMX(3, J,2)+Td)> (TCd) -PK 3,J)-TB(U TAU(3, J) ) 
£LTk‘^X(-+, J,2) =CLTRMX(4, J,2J+T(1 )'‘(TC(2)-"PI(3fJ)-TB(2) TAU(3,J)J 
GOT 080 

75 ?LTR.-)X( 1, J,2)-ELTRMXd, J,2)+T(i) ^ (-TB (1 J-* P I ( 3 , J ) -*-TC( I T AU { 3 , J J ) 

cLTRrMX(2t Ji2)=ELTRMXI 2,J, 2}+Td) • (~TB( 2) -'P I ( 3,J) + TC(2)'TAU(3,Jd 
cLTRMXI 3» J,2) =£ LTRMX ( 3 , J , 2 ) +T ( ( -TC d )“P I (3,J)+TB(1)'TAU(3,J)) 

lLTRMX( 4, J,2)=ELTRMX(4, J,2)+T(l)‘ ( -TC ( 2 P I ( 3, J I+TB ( 2 f T AU( 3 , J ) ) 
80 CONTINUE 

IF( T(4I ,LT. l.OD-14 ) GO TO 130 
N — N i 
DO 93 J = 1, JX 
PUl, J) = PI(2, J) 

PH2 , JJ = PI (3, JJ 
TAJ( 1,0) = TAU( 2, J) 

TAU(2, J) = TAU(3, J) 

90 CONTINUE 

FNAP = FNA 
FNBP = FN8 

IF (N .LE. NNX2) GO TO 65 
WRIT£(6, 9) 
wRN = .TRUE. 

RETURN 

100 DrU20J=l,JX 
C0123K=1,2 
D(U15I = 1,4 
T( I )=ELTRPXd,J,K) 

115 CONTINUE 
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£LTRMX(2, J,K) = T(l)»=2 + T(2)f-*2 
SLTRMXU,J,K) == Jl3)'*f2 + T(4)^‘#2 
ELTRMXO, J,K) = T{l)-rT{3) + T(2)^T(4) 
ELTRMXC4, J,K) = T(2J^T{3) - T(4)-^T{1) 
120 CONTINUE 

r( 1 ) = 2.0D0 '' RX-^-' 2 

QEXT = UEXT ' T(l) 

OSCAT = OSCAT •> T{1) 

CTbROS = 2.000 CTBRQS T(l) 

RETURN 

END 
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DIST 


FUNCTION DIST{RAC,A) 

REALMS A(20) ,RAC 
REALMS DlST,BfC 
B=-A( 3) 

C=RAO^*A{4) 

C=B*C 

DIST=A(U*(RAD**A(2) >^OEXP(C) 

RETURN 

END 


ORIGINAL PAGE IS 
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DIST2 


FUNCTION OIST2(R^DfB) 
REAL*8 B(20) ,RAC 
REALMS DIST2,A 
A=-(B|2)+1) 

OIST2*B{ l)*B(2)*(RAD*>»A ) 

RETURN 

END 
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APPENDIX D 


PROGRAM LISTING 
THEORETICAL SIZE 


FOR CURFIT ROUTINE USED TO FIT THE 
DISTRIBUTIONS TO THE EMPERICAL DATA 
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MAIN 


SUBROUTINE CUftFIT 

MAKES A LEAST SQUARES FIT TO A NCN-LINEAR FUNCTION 

OESCRIPTICN OF PARAMETERS 
X -ARRAY OF INC. VARIABLE CATA POINTS 

Y -ARRAY OF DSP. VARIABLE DATA POINTS 

SIGMAY -ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS 
NPTS -NUMBER GF DATA PCINTS . 

NTERMS -NUMBER OF PARAMETERS 

MODE -DETERMINES XEIGHTING FOR LEAST SQUARES FIT 
+1( INSTRUMENTAL! W(I)=1./SIGMAY( I )**2 
OINO WEIGHTIN6)W{I»=l. 

-L (STATISTICAL) WtI)=l,/Y(I) 

A -ARRAY OF PARAMETERS 

DELTAA -ARRAY OF INCREMENTS FOR PARAMETERS 
FLAMDA -PROPORTION CF GRADIENT SEARCH INCLUDED 
YFIT -array of CALCULATED VALUES OF Y 
CHISQR -REDUCED CHI SQUARE FOR FIT 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
FUNCTN(Xf I,A) 

EVALUATES THE FITTING FUNCTICN FOR THE ITH TERM 
SSP ROUTINE DSINV 

INVERTS CURVATURE MATRIX 

COMMENTS 

DATA FORMAT 

NPTStNTERMSf M0CE(3I5) 

xm,Y( I ),(SIGMAY(I) ),12(3)EL2.6) 

DIMENSION X(IOO) tY( lOC ) , S IGMAY { I 00) ,A(20 ) ,DE LTAA ( 20) ,SIGMAA(20) , 
IYFIT(100),YFITI (100) 

LOGICAL GBAD,CUR»GRID 
21 FGRMAT(3L5) 

R-EAD(5,21) GRAD, CLR, GRID 
READ(5,1) NPTS, NTERMS, MODE 

1 F0RMAT(3I5) 

IF (MCDE) 2,2,4 

2 REA0(5,3) (X(I),Y(I), 1=1, NPTS) 

3 F0RM4T(2E12.6) 

GO TC 6 

4 PEA0(5,5) ( X{ I),Y(I ),SIGHAY(I) ,1=1, NPTS) 

5 F0RMAT(3E12.6) nornrxT* 

6 READ(5,7) ( A C J) , CELTAAI J ) , J = 1 ,NT£RHS ) AL PAGE IS 

7 FORMAT(2E12.6) POOR OUAr.FTV' 

ISUM=0 ^ 

CHISC1=1.0 

14 FLAMDA=.C01 
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MAIN 


IF(CUR) GC TO 22 
IF(GRID) GO TO 23 

CALL 6RADLS(X,Y,SIGMAY,NPTSf MERMS,MCDEt A,DELTAA» 

1YFIT,CH1SQR> 

GO TO 24 

22 CALL CURFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A»OELTAA,SIGMAA,FLAMDA, 
lYFITtCHI SQR) 

GO TO 24 

23 CALL GRIDLS{X,Y ,SIGMAV,NFTS,NTERMS,MCDE,A,DeLTAA, 

1SI6MAA, YFIT.CHISQR) 

GO TO 24 

24 PRINT 8t (A(J) , J=1^NTERMSJ 

8 F0RMAT{» ',E12.6) 

PRINT RfCHISCR 

9 FORMAT! • • , ’CHISQR=« , IX, E12.6,/ J 
IF ICHISQl-CHISQR} 12,13,12 

12 CHISG1=CHISQR 
ISUH=ISUMtl 

IF ( ISUM-10) 14,13, 13 

13 DO 11 I=1,NPTS 

11 YFITHI > = 1./YFIT(I) 

PRINT 10 

10 FORMAT! • •, 13X, 'INO.VAR.* ,12X, • DEP.VAR. • , IIX , 'INV.DEP .VAR » 
PRINT 15,!Xm,YFIT(n,YFITI(I) ,I=1,NPTS) 

15 FORMAT! • ', 10X,£12.6,8X,E12.6,8X,E12.6) 

STOP 

END 
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CURFIT 


SUBROUTINE CUR F IT ( X , Y ,S IGKA V, N PTS,NTERMS ,M00 E ,A 
1,DELTAA,S1GMAA,FLAMCA,YFIT,CHSQR) 

DOUBLE PRECISION ARRAY 

DIMENSION X( 100),Y(100) , $ IGHAY ( 100 ) , A ( 2 0) ,DELTAA<20) , SIGMAAI 20 ) , 
lYFIT(lOO) ,WEIGHT< 100) tALPHA(20,20) , BETA { 20 ) , DER I V (20 ) , ARRAY ( 20 , 
120),B(20) 

11 NFREE==NPTS-NTERMS 
IF (NFREE) 13,l3-»20 
13 CHISQR=0. 

GO TC 110 

C EVALUATE WEIGHTS 

20 DO 30 I = l,NPTS 

21 IF (MCDE) 22,27,29 

22 IF (Y(D) 25,27,23 

23 WEIGHT! I ) = i,/Y( I ) 

GO TO 30 

25 WEIGHT! I )=l./(-Y( I) ) 

GC TC 30 
27 WEIGHT! I )=1. 

GO TO 30 

29 WEIGHT! I)=1./SIGNAY{I)*#2 

30 CONTINUE 

C EVALUATE ALPHA AND BETA MATRICES 

31 DO 34 J=1,NTERMS 
BETA! J»=0. 

DO 34 K=1,J 
34 ALPHA(J,K)=0. 

41 DO 50 I=1,NPTS 

CALL FDER1V(X,I ,A,DELTAA ,NT ERMS , DERI V ) 

00 46 J=1,NTERMS 

BETA(J)=BETA( J)+WEIGHT(I )*!Y(I )-FUNCTN(X, I, A ) )*DERIV( J ) 

DO 46 K=1,J 

46 ALPHAI J,K)=ALPHA{J,K)+W£IGHT( D^'DERIVfJ )*DERIV!K) 

50 CONTINUE 

51 DO 53 J=1,NTERMS 
DO 53 K=1,J 

53 ALPHA(K,J)=ALPHA!J,K) 

C EVALUATE CHISQR AT STARTING POINT 

61 DO 62 I=1,NPTS 

62 YFIT(I)=FUNCTN(X,1, A) 

63 CHI SQ1=FCHISQ( Y, SIGMA Y,NPTS, NFREE, MODE, YFIT) 

C INVFRt-xCORVAJURc MATRIX TO FIND NEW PARAMETERS 

71 DO/74 J = 1,.NTERM‘S 

72 00 73’ K=1..,NJERMS 

73 ARRAY! J,K)= ALPHA !J,K) /SORT! ALPHA! JtJJ ♦ALPHA{K,K) ) 

74 ARRAY! J,J)=1 .+FLAMDA 

80 CALL MATINV(ARRAY,NTERMS,1) 

81 00 84 J=1,NTERHS 
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BC J)-A( J) 

00 84 K=l,NTERMS 

84 B(J) = e( J)+8ETAt K)«ARRAY( J,K)/SQRT(ALPHA( J,J)'»'ALPhA(K,K) ) 
C IF CHI SQUARE I NCRE ASEO, I NCR EAS E FLAMOA 

91 DO 92 1=^1, NPTS 

92 YFIT{T)^FUNCTN(X,I,B) 

93 CHISQR=FCHISQ(YtSIGMAY,NPTS,NFREE,MOOe, YFITI 
IF (CHISCl-CHISQR) 95,101,101 

95 FLAMDA=10.*FLAM0A 
GO TO 71 

101 00 1C3 J*1,NTERHS 
103 A(JI=B(J) 

FLAMDA=FLAM0A/10. 

110 RETURN 
END 


ORIGINAL PAGE IS 
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FCERIV 


SUBROUTINE FDERI V(X, I »A,CELTAA ,NTERMS,DERIV) 
DIMENSION xaoO) ,A(20),DELTAA(20),DERIV<20) 
11 DO 18 J’slfNTERMS 
AJ=A{ J) 

DELTA=OELTAA( J) 

AU )=AJ+DELTA 
YFIT=FUNCTN{X,I , A) 

A( J)=AJ-DELTA 

DERIV(J)=(YFIT-FUNCTN(X, I,A))/(2.*DELTA) 

18 A(J) = AJ 
RETURN 
END 
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HATINV 


SUBROUTINE MATINV (ARRAY , NTE RMS i MCOOE ) 
DOUBLE PRECISION ARRAY, B 
OIMENSICN ARRAY(20,20) ,B(210) 

DO 1 I=1,NTERMS 
DO I J=L,NTERMS 

CALL LCC( I, J,I4,NTERMS,NTERMS,MCCDE) 

1 8(IJ)=ARRAY(I,J> 

EPS=l.0E-16 

CALL DSINV( BiNTERf'SfEPS, lER) 

IF (lERl 2,4,3 

2 PRINT 10 

10 FCRMAT(* *,*N0 RESULT',/) 

GO TO 4 

3 PRINT II 

11 FORMAT! • ', 'WARNING' ,/) 

4 DO 5 I=l,NTERMS 
DO 5 J=l,NTERMS 

CALL LOC( I, J,I J,NTERMS,NTERPS,MCCDE) 

5 ARRAY! I , J)=B(IJ ) 

RETURN 

END 


original 
OP POOR GiOAI"* 
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FCHISQ 


FUNCTICN FCHISQ (Y tSIGf'AY, NPTS, NFftEE, MODE, YFI T) 
DIMENSION Y( 100) ,SIGMAY( 100 l,YF IT( 100 ) 

SUM=0« 

00 5 I=1,NPTS 
IFCMODE) 1,2,3 

1 W=l./Y(I) 

GO TO 4 

2 W=l. 

GO TO 4 

3 W=l./(SIGMAY(I)**2) 

4 SUM=( Y(n-YFlf{I))*{Y(I)-YFIT( I) )+W 

5 CONTINUE 
FCHISQ=SUM/NFREE 
RETURN 

END 


.V' ’ ' 

''' ~f 


* * 
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SUBROUTINE: OSINV 


PURPOSE 

INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX 
USAGE 

CALL DSINV{A,N,EPS, lER) 

OcSCRI PTION of PARAMETERS 

A - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN 

SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT 
MATRIX. 

ON RETURN A CONTAINS THE RESULTANT UPPER 
TRIANGULAR MATRIX IN DOUBLE PRECISION. 

N - THP number OF ROWS (COLUMNS) IN GIVEN MATRIX. 

EPS - SINGLE PRECISION INPUT CONSTANT WHICH IS USED 

AS RELATIVE TOLERANCE FOR TEST ON LOSS OF 
SIGNIFICANCE. 

lER - RESULTING ERROR PARAMETER CUDEO AS FOLLOWS 
IER=J - NO ERROR 

1ER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME- 
TER N OR BECAUSE SOME RADICANO IS NON- 
POSITIVE (MATRIX A IS NOT POSITIVE 
DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI- 
FICANCE) 

IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFI- 
CANCE. THE RADICAND FORMED AT FACTORIZA- 
TION STEP K+1 WAS STILL POSITIVE BUT NO 
LONGER GREATER THAN ABS ( EPS ' A (K+ 1,K+ 1 ) ) . 

REMARKS 

THE UPPER triangular PART OF GIVEN MATRIX IS ASSUMED TO BE 
STORED COLUMNWISE IN N*{N+l)/2 SUCCESSIVE STORAGE LOCATIONS. 
IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU- 
LAR MATRIX IS STORED COLUMNWISE TOO, 

THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL 
CALCULATED RADICANDS ARE POSITIVE. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REuUIREO 
DMFSD 

METHOD 

SOLUTION IS DONE USING FACTORIZATION BY SUBROUTINE DMFSD. 


131 


ORIGINAL PAGE IS 
OF POOR QUALITT 



SUBROUTINE OSINV{4»NtEPS, lER) 


DIMENSION A(210> 

DOUBLE PRECISION A,DlNtWQRK 

FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE DMF SD 

A - TRANSPOSE(T) * T 
CALL CMFSO(A,N,EPS, lER) 

IF(IER) 9,1,1 

INVERT UPPER TRIANGULAR MATRIX T 

PREPARE INVERSION-LOOP 
IPI V = N-« ( N+l) /2 
IND = I PIV 

INITIALIZE INVERSION-LOOP 
UO 6 I=l,N 
01 N=l.DO/A(IPIV) 

A( IPIV)=DIN 
MI N=N 
KEN0=I-1 
LANF = N-KEND 
IF(KEND) 5,5,2 
J= I NO 


INITI.ALIZE ROW-LOOP 
00 A K=i,KEND 
WORK-^D.OD 
MIN=MIN-1 
LHOR=I'PIV 
LVER^J 

START INNER LOOP 
DO 3 L=LAKF,MIN 
LVER=LVER+1 
LHOk=LHGR+L 

WORK=WORK+A(LVER)' AILHORJ 
END OF INNER LOOP 

AU»=-WORK>DIN 

J=J-MIN 

END OF ROW-LOOP 

IPIV=IPI V-MIN 
IND=IND-l 

END OF INVERSION-LOOP 



on oo oo oo oooo 


OSINV 


CALCULATE INVERSECA) BY MEANS OF INVERSE(T1 
INVERSE(A) = INVERSEd) TRANSPOSE II NVERSE ( T ) } 
INITIALIZE MULTIPLICATION-LOOP 
DO 8 I=1»N 
IPIV^IPIV+I 
J=IPIV 

INITIALIZE ROW-LOOP 
DO a K=I ,N 
WQRK=0.00 
LH0R = J 

START INNER LOOP 
DO 7 L=K,N 
LVeR=LHDR+K-I 
h□RK=WQRK+A{LHOR)■^A(LVER> 

7 LHOR=LHOR+L 

END OF INNER LOOP 

A( J )=WQRK 

8 J=J+K 

END OF ROW- AND MULTIPLICATION-LOOP 

9 RETURN 
END 




of 
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SUBROUTINE DMFSO 


c 

c 

c 

c 

r 

v-> 

c 

c 

c 

c 

c 

c 

c 

c 


c 

r 

r 

C 


C 

C 


PURPOSE 

FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX 


USAGE 

CALL DMFSD{A,N,EPS»IER) 

DESCRIPTION OF PARAMETERS 

A - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN 

SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT 
MATRIX. 

ON RETURN A CONTAINS THE RESULTANT UPPER 
TRIANGULAR MATRIX IN DOUBLE PRECISION. 

N - THF NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX. 

EPS - SINGLE PRECISION INPUT CONSTANT WHICH IS USED 

AS relative TOLERANCE FOR TEST ON LOSS OF 
SIGNIFICANCE. 

lEk - RESULTING ERROR PARAMETER CODED AS FOLLOWS 
IER=0 - NC ERROR 

IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME- 
TER N OR BECAUSE SOME RADICAND IS NON- 
POSITIVE (MATRIX A IS NUT POSITIVE 
DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI- 
FICANCE) 

IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFI- 
CANCE. THE RADICAND FORMED AT FACTORIZA- 
TION STEP K^H WAS STILL POSITIVE BUT NO 
LONGER GREATER THAN ABS( EPSfA (K+ 1,K+1 ) ) . 


REMARKS 

THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE 
STORED COLUMNWISE IN N-^(N-H)/2 SUCCESSIVE STORAGE LOCATIONS. 
IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU- 
LAR MATRIX IS STORED COLUMNWISE TOO. 

THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL 
CALCULATED RADICANDS ARE POSITIVE. 

THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE 
i , SQUARE-ROOT OF THE DETERMINANT OF THE GIVEN MATRIX. 

^ /r/ ' ' If 

'•SUBIjtOOTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE* 


METHOD 

SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY, 
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THE GIVEN MATRIX IS REPRESENTED AS PRODUCT OF TWO TRIANGULAR 
MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF 
THE RETURNED RIGHT HAND FACTOR. 


SUriROUT INE OMFSD(A,N, CPS, lER) 


DIMENSION A(210> 

DOURL6 PRECISION OPIV,DSUM,A 

TEST ON WRUNG INPUT PARAMETER N 
IF(N-l) 12,1,L 

1 IER=0 

INITIALIZE DIAGONAL-LOOP 
API V=0 

DO 11 K=1,N 
KPIV=KPIV+K 
IND=KPIV 
LEND=K-1 

CALCULATE TOLERANCE 
TOL=A6S(EPS SNGLt A(KPIV) n 

START FACTORIZATION-LOOP OVER K-TH ROW 
DO 11 I=K,N 
DSUM=0.00 
IFILENDJ 2,4,2 

START INNER LOOP 

2 DO 3 L=1,L£ND 
LANF=KPIV-L 
LINO=IND-L 

3 DSUM=OSUM*-A(LANF) rAIL IND) 

C END OF INNER LOOP 

C 

C TRANSFORM ELEMENT A(INO) 

^ QSUM=A( IND)-DSUM 
IF(I-K) 10,5,10 
C 

C TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SI6NIFICANCF 

5 IF(SNGL{OSUM)-TOL) 6,6,9 

6 iriDSUM) 12,12,7 

7 IF(IER) 8,8,9 

8 IER=K-1 
r 
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DMFSD 


C COMPUTE PIVOT ELEMENT 

9 DPIV=DSQRT{DSUM) 
A{KPIV)=DPIV 
DP I V=l. DO/ DP IV 
GO TO 11 
C 

C CALCULATE TERMS IN ROW 

10 AUND) = OSUM*'DPIV 

11 IND=IND+I 
END OF DIAGONAL-LOOP 

RETURN 

12 IER=-1 
RETURN 
END 
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APPENDIX E 


RELATIONSHIP BETWEEN EXTINCTION, SCATTERING, AND 
ABSORPTION COEFFICIENTS AND THE MLE PARAMETERS 
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The extinction (a), scattering (s), and absorption (a) coeffi- 
cients for suspended particulates can be calculated using the Mie 
formalism. Using the Mie parameters, an(x,m) and bn(x,m) of 
equations (3-7) the extinction coefficient is given by: 


a. = 2^ I (2n+l) ^Re(an(x,m)) + Re(bji(x,m))^ n(r)dr (E-1) 


2it 


n=l 


where n(r) is the particle size distribution function and x = 2-rrr/X. 
The expression for the scattering coefficient is: 


2 /' <» 

s = X / (2n+l) 


2ir 


E 

n=l 


{|a„ 


(x,m) ^ + jbn(x,m)|^ J n(r)dr (E-2) 


The absorption coefficient is the difference between a and s, thus 


a 


X^y~^ (2n+l) -j^ReCanCx.m)) + Re(bn(x,m)) - 

|a^(x,m)j^ - |bn(x,in) 1 ^ 1 n(r) dr. 


(E-3) 


The A7a-lues for a, s, and a used xn the Monte Carlo routine were not 

’ ^ ^ f i » it t * 

calculated'* in I this way because the values explicitly depend on the 

^ /i , 

concentration through n(r)* Instead a, s, and a were chosen to 

correspond to physically observed values. 

The absorption coefficient depends on the xmaginary part of 

the index of refraction, but in a non- trivial way. If Im(m) == 0 

(25) 

then it can be shown ^ that 
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(E-4) 


j a^(x,m) - %( == 3a 

2 

jbn(x,in) - %j = 3a . 

Expanding equation (E-4) leads to 

r 2 2 

[Re(an(x,m))3 - Re(an(x,m)) + [lm(an(x,m))] + 3a = ^ 


2 r 2 

Re(an(x,in)) = [Re(an(x,m))J + [lm(an(x,m) )] 

= an(x,m) ^ (E-5) 

with a similar result holding for bn(x,m). Using these results in 
equation (E-3) leads to a=o. Thus if the imaginary part of the 
index of refraction is zero the absorption coefficient is also zero. 
If Ira(m) ^ 0 then^^^^ 
a„(x,m) - <k 

jbn(x,m) - < k 

Which, after expansion, leads to 

I 1 2 

Re(a„(x,m)) > |a (x,m)| 

I 1 2 

Re(bn(x,m)) > |b (x,m)| , 

so that, by equation (E-3), a>o for a non-zero imaginary component 
in the index of refraction. 
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